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CHAPTER 

ONE 

Introduction  and  Systems  Approach 


In  recent  years,  the  size  of  satellites’  structures  that  are  being  developed  has 
increased  dramatically.  Operational  power  requirement  growth  has  outpaced 
the  ability  of  standard  solar  panels.  Correspondingly,  new  design  problems  have 
arisen.  In  order  to  handle  these  technological  difficulties,  a  systematic  analysis 
of  the  space  structures  is  required.  This  study  incorporates  a  cross-disciplinary 
approach  that  is  required  in  order  to  analyze  a  technological  requirement,  namely 
on-orbit  structural  identification. 

In  the  following  sections,  the  basis  of  the  need  for  larger  satellites  and  so¬ 
lar  panels  is  explained.  The  power  requirement  solution  may  follow  different 
technical  paths.  A  very  simple  trade-off  analysis  leaves  the  designer  with  solar 
panels  as  the  preferred  power  source.  The  power  requirement  thus  drives  the 
requirement  for  larger  sized  panels.  These  larger  panels  cause  structural  prob¬ 
lems  which  require  monitoring  and  analysis.  The  identification  and  monitoring 
requirement  leads  to  the  topic  of  this  thesis,  which  is  ground  based  laser  mea¬ 
surement  of  satellite  motion  and  the  parameter  identification  of  the  structural 
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modes. 


A  top-down  diagram  of  the  technological  flow  that  drives  this  experiment 
is  shown  in  Figure  1.1.  This  figure  diagrams  how  disperse  top-level  factors  of 
requirements  lead  to  the  need  for  large  space  systems.  These  systems  generate 
technological  concerns,  one  of  which  is  structural  analysis.  The  solutions  of  this 
concern  fan  out  to  a  variety  of  choices.  The  choice  of  remote  sensing  is  the  topic 
of  this  thesis.  The  following  sections  will  explain  in  more  detail  the  technological 
factors  that  contribute  to  the  system.  \  - 
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1.1  Technological  Forces  on  Spacecraft  Size 


There  are  two  technological  forces  driving  the  increasing  size  of  space  structures. 
The  first  is  the  increased  on-board  power  demands  for  communication,  scientific, 
and  military  payloads.  Typical  communication  satellites  require  between  1  and 
1.5  kilowatts  of  power  during  their  lifetimes.  The  new  SDI  space-borne  weapon 
systems  may  require  20  times  this  amount.  These  power  requirements  have 
surpassed  the  capacity  of  solar  panels  directly  mounted  to  the  rigid  body.  To 
satisfy  these  requirements  either  solar  panels  or  nuclear  reactors  can  be  used  (see 
Figure  1.2).  The  latter  source  is  typically  an  unacceptable  solution  because  of 
safety  and  political  reasons.  Still,  the  current  technology  for  solar  panels  can 


Figure  1.2:  Modern  Satellite  Power  Sources 


only  keep  up  with  such  demands  by  increasing  the  size  of  solar  panels.  Simply 
increasing  the  size  of  the  spacecraft,  however,  will  not  provide  enough  surface 
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area  exposed  to  solar  radiation  to  meet  the  required  power  demands. 

The  solution  has  been  to  use  large  deployable  solar  panels  or  wings.  These 
large  wings  are  light  weight  and  constructed  with  graphite  composite  material 
but  are  susceptible  to  large  flexible  motion.  This  motion  may  be  driven  by 
solar  wind,  atmospheric  drag,  or  interaction  with  the  spacecraft’s  main  body. 
Rigid  body  dynamics,  system  vibratory  motion,  and  solar  tracking  actuators 
characterize  this  latter  motion  which  is  transferred  to  the  solar  wings  at  the 
wing/bus  interface  thus  exciting  the  wings’  vibration  (See  Figure  1.3). 

The  second  technological  force  is  the  development  of  complex  scientific  ex¬ 
periments  and  payloads.  These  payload  require  large  structures  such  as  the 
space  station  or  complex  antenna  booms.  Construction  of  such  structure  will 
require  extended  light  weight  booms  and  elements.  These  elements  typically  will 
be  built  from  hinged  graphite  elements  and  deployable  truss  booms  (See  Figure 
1.4).  The  spacecraft  that  is  the  subject  of  this  study  will  fall  into  this  second 

category  of  large  space  structures. 

These  new  space  structures  utilize  light  weight  materials  which  will  have  low 
damping  factors.  Stiffness  requirements  for  structural  rigidity  contribute  to  the 
low  internal  damping.  The  combination  of  large  elements  and  low  damping  cre¬ 
ates  new  design  challenges  as  the  spacecraft  lowest  vibration  frequencies  begin  to 
extend  into  the  lower  frequencies.  Because  of  light  damping,  all  control  signals 
or  input  forces  are  required  not  to  generate  signal  with  frequency  components 
near  these  vibrational  modes.  Any  dynamic  motion  that  generates  forces  near 
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these  modes  may  excite  large  vibrational  motion.  This  resonant  motion  may 
present  an  operational  nuisance,  but  in  some  cases,  it  could  become  disastrous. 
Unfortunately,  as  the  scale  of  deployed  structures  get  larger,  the  lower  modal 
frequencies  may  begin  to  seriously  limit  the  actuator  operating  range.  Addi¬ 
tionally,  many  actuators  provide  incremental  pulses  (such  as  stepper  motors 
and  thrusters)  and  the  spectrum  of  such  pulses  are  broad  and  will  excite  the 
structural  modes. 


Figure  1.3:  Modern  satellite  high  power  solar  panels 
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Figure  1.4:  Modern  space  structures 


1.2  The  Need  for  Structural  Identification 

While  structural  analysis  can  predict  modal  frequencies  with  some  degree  of  con¬ 
fidence,  the  structure  configuration  over  a  space  structures  life-cycle  will  not  be 
constant.  Normal  reconfiguration  requirements,  mechanism  decay  and  failure, 
and  changes  in  structural  properties  will  create  complex  structural  configura¬ 
tions  that  may  be  increasingly  hard  to  model,  especially  as  unknown  elements 
creep  into  the  model  formulation.  Therefore,  there  exists  a  need  to  indepen¬ 
dently  verify  and  identify  structural  behavior.  A  structural  identification  of  the 
spacecraft  in  its  “true”  environment  will  provide  confidence  in  analytical  models 
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and  hopefully  point  to  areas  where  more  research  is  needed.  Identification  can 
be  used  to  monitor  the  spacecraft  over  its  life  and  analyze  trends  in  performance. 
Finally,  in-flight  modeling  is  a  major  step  toward  the  goal  of  real-time  adaptive 
attitude  contol  and  stabilization. 

Interest  in  structural  identification  has  recently  expanded.  Tests  have  been 
done  by  the  US  Air  Force  that  tested  structural  modes  on  a  boom  while  in  a 
zero-gravity  environment.  This  was  done  in  the  hold  of  an  airplane  while  in  a 
parabolic  dive.  Currently,  there  are  plans  to  perform  a  similar  experiment  m 
the  hold  of  the  space  shuttle.  The  experiment  of  this  thesis  has  the  advantage 
of  testing  flight  hardware  using  simple  low  cost  hardware  and  accomplishing 

similar  goals. 

1.3  Thesis  Overview 

This  thesis  accomplishes  parameter  identification  on  an  m-orbit  satellite.  In 
Chapter  2,  the  satellite  on  which  the  structural  identification  will  be  performed 
is  described.  In  Chapter  3,  the  structural  modeling  techniques  used  to  model 
the  spacecraft  are  presented.  In  Chapter  4,  various  algorithmic  approaches  are 
described  that  handle  identification  when  dealing  with  noise  corrupted  data. 
In  Chapter  5  the  actual  experiment  is  described  and  the  results  are  presented. 
These  results  are  fascinating  in  that  using  only  small  observation  windows,  the 
first  three  of  four  structural  modes  were  detected. 
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CHAPTER 

TWO 

Spacecraft  Configuration  and  Experiment  Description 


The  experiment  was  designed  for  the  Low  Atmospheric  Control  Experiment 
or  LACE  spacecraft.  It  is  a  low  earth  orbit  satellite  that  was  launched  into  a 
circular  orbit  with  an  altitude  of  approximately  550  kilometers  and  an  inclination 
of  43  degrees.  The  LACE  spacecraft  structure  is  composed  of  a  central  rigid  body 
or  bus  with  three  deployable  booms  (See  Figure  2.1).  The  bus  carries  mission 
sensors  and  experiments,  all  supporting  telemetry/command  modules,  attitude 

and  control  subsystem,  and  the  solar  panels. 

Each  of  the  deployable  booms  has  a  different  mission  function.  The  first 
boom  is  the  gravity  gradient  boom  and  is  oriented  directly  away  from  the  earth. 
It  has  an  electro-magnetic  energy  dissipating  unit  located  at  the  tip.  The  energy 
dissipation  unit  is  part  of  a  passive  attitude  stabilization  system  that  dumps 
destabilizing  dynamic  energy  using  the’ earth’s  magnetic  fields. 

The  second  boom  is  the  retro-reflector  boom  and  is  deployed  in  the  direc¬ 
tion  of  the  spacecraft  velocity  vector.  The  retro-reflector  boom  tip  has  a  laser 
reflector  unit  mounted  on  it.  These  reflectors  are  part  of  the  primary  spacecraft 
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mission.  Amongst  these  reflectors  is  a  germanium  reflector  (approximately  1 
inch  diameter)  that  is  dedicated  to  this' structural  dynamics  experiment. 

The  third  boom  is  the  balance  boom  and  is  oriented  180  degrees  from  the 
retro-reflector  boom.  The  balance  boom  has  a  strictly  passive  role  of  counter¬ 
acting  the  rigid  body  dynamics  due  to  the  retro-reflector  boom.  The  rigid  body 
and  the  balance  boom  also  have  germanium  reflectors  mounted  upon  them. 

The  spacecraft  booms  are  both  deployable  and  retractable.  The  deployed 
length  varies  from  0  feet  up  to  150  feet.  The  booms  are  constructed  of  light 
weight  composite  material  with  continuous  longerons  and  stiff  cross  elements 
(See  Figure  2.2).  Since  the  boom  lengths  are  variable,  the  system  vibration 
modes  are  variable  as  well  and  are  a  function  of  the  deployed  length  .  The  un¬ 
deployed  portion  of  the  booms  remains  elastically  coiled  in  the  deploying  canister 
mounted  to  the  main  spacecraft  body.  Additionally,  the  deploying  canister  has 
an  elastic  compliance.  This  compliance  is  incorporated  in  the  vibration  analysis, 
see  Table  2.1  (also  see  Chapter  3). 


density 

unit 

length 

kg/m 

bending 

stiffness 

(El) 

N-m2 

torsional 

stiffness 

(GJ) 

N-m2 

Rotating 
Inertia  I0 
unit /length 
kg-m 

boom  canister 
compliance 

N /radian 

0.291 

«  1.55  x  104 

631  1-33  x  10 

1.695  x  104 

Table  2.1:  Boom  Structural  Properties 


The  rigid  central  body  has  a  mass  of  1177  kilograms  with  moments  and 
products  of  inertia  listed  in  Table  2.2.  These  properties  were  experimentally 
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measured.  The  tip  masses  on 


each  boom  are  listed  in  Table  2.3. 


Inertia 

kg  -  m2 

I XX 

144S.7 

Ly 

1426.4 

Lz 

1026.2 

Ly 

3.61 

Lx 

19.98 

Iyz 

14.S6 

Table  2.2:  LACE  Bus  Inertial  Properties 


boom 

kg 

!  Gravity 

90.7 

Retro-reflector 

15.9 

Balance 

15.9 

Table  2.3:  LACE  Boom  tip  masses 
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Figure  2.1:  Spacecraft  Configuration 
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Figure  2.2:  Deployable  boom  construction 
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2.1  Experiment  Description 


The  experiment  is  to  measure  linear  vibration  of  the  complete  structure  and 
calculate  the  system  modal  frequencies  and  damping  factors.  The  vibrations 
will  be  excited  by  the  deployment  and  retraction  of  the  lead  boom.  This  type  of 
excitation  is  due  to  the  nonlinear  coupling  of  the  elastic  strains  of  the  booms  in 
the  canister.  Additionally,  the  lead  tip  mass  is  offset  from  the  axis  of  the  boom. 
This  offset  provides  a  moment  at  the  boom  tip  which  will  excite  vibration.  While 
the  excitation  forces  involve  nonlinear  factors,  the  resulting  vibration  will  end 
up  being  within  the  linear  range  (See  chapter  3  on  beam  modeling). 

Sampling  the  vibrations  will  be  done  remotely.  Since  the  LACE  spacecraft 
does  not  have  telemetered  data  that  directly  indicates  vibratory  motion  (  such 
as  strain  gauges,  gyros),  the  measurements  must  be  taken  indirectly  or  using 
off-board  devices.  Originally,  the  experiment  had  been  designed  to  use  two 
techniques  to  acquire  vibration  data. 

1.  The  first  technique  uses  a  ground  based  laser  to  illuminate  the  germanium 
mirrors  on  the  satellite  (See  section  1.1).  The  reflected  laser  signal  will 
have  a  Doppler  shift  due  to  the  relative  motion  of  the  satellite  with  respect 
to  the  laser  source.  The  detected  Doppler  shift  from  each  of  the  illuminated 
mirrors  is  then  used  to  calculate  this  relative  motion.  (See  Figure  2.3) 

2.  The  second  technique  uses  a  spacecraft  mounted  optical  sensor  to  detect 
vibration  motion.  The  spacecraft  sensor  used  to  track  stars  or  ground 
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based  lasers.  This  sensor  can  track  within  precise  toleiances  (  ~  millira- 
dians  ).  Using  this  sensor,  the  vibrations  that  influence  the  central  body 
can  be  detected.  Thus  by  calculating  several  coordinate  transformations, 
the  motion  detected  by  the  sensor  can  be  projected  onto  the  expected 
vibration  axis. 

Although  considerable  effort  was  put  into  the  second  approach,  it  was  eliminated 
eventually  from  the  experimental  technique  due  to  operational  considerations. 
Specifically,  there  exists  an  operational  constraint  against  optical  tracking  while 
the  spacecraft  was  being  illuminated  by  an  800  watt  infrared  laser.  Considering 
that  the  reflected  laser  return  provides  greater  resolution,  the  sensor  tracking 
approach  had  to  be  eliminated. 

2.1.1  Doppler  Shifted  Laser  Return 

The  laser  source  and  detection  facility  used  is  located  at  Lincoln  Laboratories 
in  Massachusetts.  This  facility  was  contracted  to  perform  the  actual  laser  mea¬ 
surements.  The  collected  data  was  then  sent  to  the  Naval  Research  Laboratory 
where  all  the  analysis  was  performed. 

The  data  collection  is  done  by  pulsing  the  spacecraft  with  an  infrared  laser 
source  and  measuring  the  Doppler  shifted  return.  The  laser  pulsing  used  a  duty 
cycle  of  62  Hz.  This  pulsing  rate  is  higher  than  the  modes  targeted  for  detection 
(Ref.  Chapter  3),  but  the  oversampling  provides  a  more  reliable  data  set  from 
which  noise  is  more  easily  eliminated.  Each  pulse  is  sampled  40S0  times  in-phase 
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and  in-quadrature  (1.25  Mhz  sample  rate).  The  frequency  components  of  the 
Fourier  transform  of  the  complex  data  set  is  analyzed  to  detect  Doppler  shifts  in 
the  return  spectrum  (See  Figure  2.4).  The  Doppler  shifts  of  each  return  source 
is  calculated  from  the  sampled  power  spectrum.  The  relative  Doppler  shift  of 
each  boom  compared  to  the  central  bus  is  used  to  determine  vibratory  motion 

(See  Figure  2.5). 

In  addition  to  the  vibration  motion,  the  Doppler  shift  will  include  the  rigid 
body  motion  of  the  satellite’s  orbital  motion.  This  rigid  body  motion  becomes 
significant  on  the  line  of  sight  vector  from  an  earth  base  observation  reference. 
(See  Figure  2.6).  Considering  that  observation  windows  are  time  limited,  elimi¬ 
nation  of  this  rigid  body  component  will  improve  the  potential  resolution  of  the 
lowest  vibration  frequency  modes  [15]  (Reference  Chapter  3  for  modal  frequency 
analysis). 
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shill  In  kHz 


Delta  v  in  kiloHertz:  predicted  orbit  of  14  July,  0258 


Figure  2.6:  Observed  Rigid  Body  motion  Day  258 


CHAPTER 

THREE 

Modeling  the  Spacecraft  Structural  Properties 


The  modal  analysis  can  be  performed  in  two  ways.  The  first  approach  is  finite 
element  modeling  (FEM),  and  the  second  technique  is  a  continuum  approach 
through  a  partial  differential  equation  model  (PDE).  The  LACE  geometry  mod¬ 
eling  lends  itself  to  straightforward  PDE  solution.  Once  formulated,  the  solution 
of  the  PDE  provides  a  quick  and  accurate  solution.  FEM  also  provides  conve¬ 
nient  three  dimensional  modeling,  but  its  results  will  converge  to  PDE  solutions 
only  as  the  number  of  nodes  becomes  large.  The  advantage  of  FEM  analysis 
is  that  software  packages  are  readily  available  and  easy  to  use.  For  this  work, 
PDE  solutions  were  used  and  FEM  solutions  were  used  as  a  comparison. 

The  PDE  model  provides  some  distinct  advantages  to  the  FEM  approach. 
One,  the  number  of  parameters  used  in  the  model  is  greatly  reduced.  For  the 
LACE  experiment,  the  number  of  unknown  structural  properties  are  the  boom 
stiffness  and  torsional  stiffness,  boom  canister  compliance,  and  boom  internal 
damping  factors  (assuming  all  booms  are  identical).  Additionally,  any  struc¬ 
tural  irregularities  are  lost  within  the  homogenized  model.  The  benefit  of  this 


19 


homogenization  process  was  theoretically  explained  by  Blankenship  [10].  Fi¬ 
nally,  continuous  models  can  be  used  directly  in  control  analysis.  The  large 
order  FEM  models  must  be  approximated  using  model  truncation  techniques 
before  control  system  analysis  can  be  performed  (  [18],  [14]). 

3.1  Euler  Beam  Model 

The  basic  structural  model  used  was  the  Euler  beam.  The  Euler  beam  is  the 
traditional  choice  because  it  is  a  linear  model.'  More  complicated  models  may 
provide  more  “accuracy”,  however,  they  introduce  nonlinear  components.  These 
nonlinear  components  obscure  the  definition  of  modal  frequency  and  damping. 
Additionally,  the  degree  of  deflection  expected  in  the  experiment  are  within  the 
engineering  bounds  of  linearity. 

The  Euler  beam  basic  assumptions  are  : 

1.  The  cross-sectional  dimensions  are  <<  than  the  length  of  the  beam.  This 
is  the  slender  beam  condition.  For  the  booms  of  the  LACE  spacecraft, 
this  condition  is  easily  satisfied. 

2.  A  plane  cross-section  drawn  normal  to  the  beam  axis  remains  planar  and 
normal  to  the  beam  axis  after  bending. 

3.  Shear  displacement  is  very  small  compared  to  bending  displacement. 

4.  The  amount  of  rotation  of  any  particular  cross-section  of  the  beam  is  << 
than  the  amount  of  deflection  of  the  cross-section. 
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While  these  assumptions  may  be  severe,  they  are  essentially  obeyed  with  the 
small  excitation  due  to  boom  deployments.  The  booms  are  assumed  to  be 
cylindrical  beams.  This  assumption  was  used  in  the  FEM  analysis  also.  Thus, 
the  PDE  model  describing  the  flexible  motion  of  the  beam  through  any  axis  is 
described  by: 


d2w(z,t)  r^Td‘lw(z,t)  _ 
p—W—  +  EI—dF~- 


(3.1.1) 


This  equation  assumes  the  material  stiffness  El  and  mass  per  unit  length  p  are 
uniform  along  the  beam.  The  values  for  El  and  p  used  in  this  equation  are  the 
same  as  those  in  Table  2.1.  The  solution  is  calculated  by  using  the  Laplace 
transform  on  the  variable  time  in  the  equation.  The  equation  (  3.1.1)  becomes: 


d4w(z,s)  s2p^/  v 

dz 4  ei  v  ' 


(3.1.2) 


The  solution  to  (  3.1.2)  is: 


w(z,  s) 
A 


A\CO$Xz  -f-  A2sinXz  -f-  A^coshXz  +  A^sinhXz 


(— )?■ 
yEI} 


(3.1.3) 

(3.1.4) 


Given  (  3.1.4)  and  boundary  conditions  at  each  end  of  the  beam,  the  co¬ 
efficients  Ai  i  =  1,  *  -  - ,  4  are  solved.  The  boundary  conditions  considered  for 
this  experiment  are  tip  masses  and  the  boom  canisters  represented  as  spring 

elements  at  the  bus  (See  Figure  3.1). 

Likewise,  the  torsional  partial  differential  equation  is  described  by. 


.  d29(z,t) 

0  <9f2 


+  GJ 


dz 2 


(3.1.5) 
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This  equation  assumes  the  torsional  stiffness  GJ  and  rotating  mass  inertia  per 
unit  lenght  I0  are  uniform  along  the  beam.  The  values  for  GJ  and  I0  used  m 
this  equation  are  the  same  as  those  in  Table  2.1.  The  solution  is  calculated  by- 
using  the  Laplace  transform  on  the  variable  time  in  the  equation.  The  equation 
(  3.1.5)  becomes: 


d2d(z,s) 
dz 2 


+  “ gj Q{Z,  S)  —  0 


(3.1.6) 


The  solution  to  (  3.1.6)  is  shown  to  be 


6(z,s) 
! 3 


Aicosfiz  + 


yGj' 


A^sinfiz 


(3.1.7) 

(3.1.8) 


Using  equation(  3.1.8)  and  the  boom’s  boundary  conditions,  the  coefficients  At- 
are  solved. 

The  complete  PDE  solution  is  calculated  be  solving  the  matrix  of  boundary 
conditions  for  all  flexible  motion  equations.  The  individual  booms  connected  at 
the  spacecraft  bus  satisfy  the  dynamic  and  boundary  conditions  of  the  bus.  The 
modal  frequencies  are  then  calculated  from  the  boundary  condition  matrix  . 
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3.2  Model  Analysis  Results 


Using  the  model  of  the  previous  section,  the  modal  properties  of  the  LACE 
spacecraft  were  determined.  The  spacecraft  boom  length  configuration  analyzed 
is  listed  in  Table  3.1.  This  was  the  configuration  of  the  spacecraft  at  the  time 
of  the  identification  experiment. 

The  pitch  modes  are  the  most  likely  detectable  modes.  For  simplification  the 
pitch  modes  were  isolated  from  the  yaw  and  roll  modes.  This  speeds  up  analysis 
considerably.  The  isolation  is  accomplished  by  assuming  the  cross  products  of 
inertia  are  zero  and  that  the  booms  are  mounted  exactly  m  line  with  the  bus 
center  of  mass.  This  assumption  is  reasonable  given  the  values  for  the  product 
of  inertia  as  shown  in  Table  2.2  are  small.  A  study  of  the  lowest  pitch  modes 
showed  that  the  pitch  modal  frequencies  were  not  sensitive  to  this  assumption. 
Of  course,  for  a  more  exact  model  all  the  vibration  planes  must  be  considered. 
The  actual  bus  dimensions  do  marginally  affect  calculated  pitch  modes.  The 
yaw  and  roll  modes  were  calculated  separately.  In  Table  3.2  the  pitch  modes 
are  listed.  Additionally,  Table  3.2  lists  the  relative  magnitude  of  displacement 
experienced  at  the  lead  tip  for  each  mode.  Relative  displacement  is  the  ratio  of 
tip  deflection  to  peak  mode  deflection.  This  ratio  gives  some  insight  to  which 
modes  are  most  likely  to  be  detected  with  the  ground  based  laser.  Displacement 
ratios  which  exhibit  large  motions  at  the  modal  frequency  are  more  likely  to  be 
detected  (See  Figure  3.2  through  3.7).  However,  it  should  be  noted  that  there 
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is  no  guarantee  any  particular  mode  will  be  excited.  This  is  due  to  the  lack 
of  control  of  the  excitation  mechanism.  It  should  be  noted  that  the  vibration 
mode  at  .7567  Hz  is  primarily  due  to  the  modeled  spring  compliance  of  the 
boom  canister.  Changing  the  boundary  condition  at  the  spacecraft  to  a  fixed 
clamp  drastically  changes  this  mode.  As  the  modal  frequencies  get  higher,  the 
displacement  ratios  get  larger.  This  suggests  that  if  the  system  higher  modes 
could  be  excited,  then  they  could  be  easily  detected.  Unfortunately,  the  low 
energy  excitation  force  used  in  this  experiment  is  not  expected  to  generate  the 
higher  modes. 


Boom 

length  (ft) 

Gravity  Gradient 

150 

Retro- Reflector 

15 

Balance 

150 

Table  3.1:  LACE  Boom  Configuration  Analyzed 


PDE 

hz 

FEM 

hz 

Displacement 

Ratio 

.01908 

.01906 

0.0242 

.1298 

.12981 

0.0063 

.2581 

.25782 

0.0284 

.3238 

.32333 

0.0565 

.7567 

.74952 

0.4669 

.8217 

.S1SS0 

0.1068 

Table  3.2:  LACE  System  Pitch  vibration  modes  below  1  Hertz 
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Figure  3.3:  Pitch  Plane  modal  shape:  .1298  Hz 
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Figure  3.5:  Pitch  Plane  modal  shape:  .3238  Hz 
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Figure  3.6:  Pitch  Plane  modal  shape:  .7567  Hz 
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CHAPTER 

FOUR 

Structural  Identification  using  Free-Decay  Observations 


The  LACE  structural  identification  problem  presents  unique  challenges.  As 
described  in  chapter  3,  the  form  of  the  excitation  or  driving  function  cannot 
be  accurately  modeled.  Thus,  the  analysis  will  be  limited  to  only  the  post¬ 
excitation  free-decay  measurements.  Additionally,  the  observable  measurement 
windows  available  for  laser  reflections  are  extremely  limited  in  time.  Given  these 
two  major  limitations,  the  identification  analysis  must  be  robust  in  the  presence 
of  significant  noise  and  must  converge  to  a  model  given  only  a  small  window  of 
measurements. 

4.1  Identification  Introduction 

The  project’s  goal  is  the  accurate  identification  of  the  LACE  spacecraft’s  system 
modal  frequencies  and,  if  possible,  damping  factors  .  Such  experiments  are  vital 
compliments  to  a  structural  analysis  computed  by  either  finite  element  models 
or  continuous  partial  differential  equation  models.  A  successful  structural  iden¬ 
tification  of  the  spacecraft  in  its  “true”  environment  will  provide  confidence  m 
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analytical  models  and  hopefully  point  to  areas  where  more  research  is  needed. 
Additionally,  in-flight  modeling  is  a  primary  component  toward  the  goal  of  real¬ 
time  adaptive  attitude  control  and  stabilization. 

As  mentioned,  the  experiment  will  rely  completely  on  free  decay  measure¬ 
ments.  Thus  limited,  the  system  modeling  will  attempt  to  duplicate  the  struc¬ 
tural  internal  dynamics  with  an  impulse  response  model.  The  identification  will 
utilize  linear  system  theory  which  provides  error  bounds  on  the  model  s  per¬ 
formance.  The  performance  criteria  that  will  be  considered  are  the  functional 
norms  I2,  X2[0,oo),  Too,  Frobenius  (  ||.||f  ),  Nuclear  (  ||.||jv  )  and  the  Hankel  ( 
||.||K  ).  The  system  identification’s  performance  is  measured  by  the  comparison 
of  the  true  structural  free-decay  and  the  linear  model’s  impulse  response.  Per¬ 
formance  in  both  the  time  domain  and  the  frequency  domain  will  be  addressed. 

Historically,  the  identification  of  large  structures  by  using  an  impulse  re¬ 
sponse  is  relatively  recent.  The  most  striking  improvement  to  structural  identi¬ 
fication  was  the  Eigensystem  Realization  Algorithm  (ERA)  suggested  by  Juang 
et  al  [13].  This  algorithm  is  theoretically  inspired  by  the  Kalman-Ho  algorithm 
[2]  but  differs  by  utilizing  the  more  numerically  robust  approach  of  singular  value 
decomposition.  Even  more  recently,  Mook  et  al  [4]  presented  the  Minimum 
Model  error  approach.  This  approach,  coupled  with  ERA,  provides  improved 

performance  over  that  of  ERA  alone  [5]. 

While  these  two  approaches  were  effective  in  parametric  identification,  they 

did  not  provide  the  theoretical  background  to  explain  the  accuracy  of  the  al- 
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gorithms  or  to  calculate  any  error  bounds  associated  with  the  identification. 
In  this  chapter  the  error  bounds  associated  with  ERA  and  ERA/MME  will  be 
shown. 

In  section  4.2  the  background  of  the  Hankel  operator  is  presented  and  the 
error  bounds  of  truncating  models  are  shown  for  balanced  truncation  and  of  opti¬ 
mal  Hankel  norm  reductions.  In  section  4.3  the  modeling  of  infinite  dimensional 
systems  will  be  addressed.  In  section  4.4  the  ERA  algorithm  is  derived  con¬ 
sidering  the  balanced  realization  approach.  In  section  4.5  the  Minimum  Model 
Error  approach  is  explained  using  Pontryagin’s  approach,  and  the  solution  to 
the  two-point  boundary  value  problem  is  solved  using  linear  algebra  as  proposed 
by  Mook  [6].  In  section  4.6  a  novel  approach  of  least-squares  modeling  with 
balanced  order  reduction  will  be  proposed.  In  section  4.7  explanation  of  opti¬ 
mal  sampling  is  discussed.  Finally,  section  4.8  an  example  will  be  explored  using 
data  generated  from  the  impulse  response  of  a  spacecraft  model. 

4.2  The  Hankel  Operator  and  Error  Bounds  of  Model  Reduction 
schemes. 

The  Hankel  operator  and  the  properties  of  its  singular  values  are  central  to  all 
model  order  reduction  schemes.  In  this  section  the  relation  of  the  Hankel  oper¬ 
ator  and  the  controllability  and  observability  grammians  is  shown.  Using  these 
properties,  a  linear  model  can  be  constructed  that  will  model  the  true  response 
within  a  given  bound.  In  the  background  development,  the  Hankel  operator  is 
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considered  for  both  discrete  and  continuous  systems.  The  motivation  for  this  is 
that  although  the  spacecraft  structural  system  is  a  distributed  continuous  infi¬ 
nite  dimensional  system,  all  the  experimental  data  is  confined  to  the  sampled 
discrete  domain.  It  was  shown  that  the  singular  values  of  a  continuous  system 
Hankel  operator  are  identical  to  its  discretized  counterpart  [11]. 

4.2.1  The  Hankel  Operator 

All  realizations  constructed  utilizing  the  properties  of  the  Hankel  operator  will  be 
modeled  using  the  standard  state  space  notation  [16].  For  finite  order  continuous 
systems  this  representation  is: 


x(t)  =  Ax(t)  +  Bu(t) 

(4.2.1) 

y(t )  =  Cx(t)  +  Du(t) 

(4.2.2) 

For  discrete  systems  this  representation  is: 

x(n  +  1)  =  Ax(n)  +  Bu{n) 
y(n)  =  Cx(n )  +  Du{n) 

Given  this  format,  the  controllability  and  observability  operators  are  then  de¬ 
fined. 

Definition  4.1  For  a  given  a  state  space  representation  ( A,B,C ),  the  control¬ 
lability  operator  maps  future  inputs  to  the  internal  states  # c  :  L2[0,cc)  -+  Cn . 
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For  continuous  systems  this  is  defined  as: 

roo 

ycu  =  -  /  exp(—Ar )Bu(r)dT 
Jo 

For  discrete  systems  \I/C  :  Z2[0,  oo)  ~ *  Cn  is  defined  as. 

oo 

ycu  =  AkBu(k) 

k=0 

Definition  4.2  For  a  given  state  space  representation  ( A,B,C )  and  internal 
state  x  6  Cn,  the  observability  operator  maps  the  internal  state  to  past  outputs 
v£o  .  Qn  ^2(_OO)0]  .  For  continuous  systems  this  is  defined  as: 

(f0i)(f)  =  Cexp(At)x 

For  discrete  systems  \l/0  :  Cn  — >  /2(— so50]  z5  defined  as. 

(90x)(k)  =  CAkx 

Note:  That  if  (A,  B)  is  controllable  then  is  surjective  and  if  (A,  C)  is  observ¬ 
able  then  $0  is  injective.  Since  these  operators  are  Hilbert  space  isomorphisms 
(ie.  they  map  from  one  Hilbert  space  to  another),  the  adjoint  operators  exist. 

Definition  4.3  For  a  given  a  state  space  representation  ( A,B,C ),  the  adjoint 
controllability  operator  maps  internal  states  to  future  inputs  :  Cn  -»  X2[0,oo) 
For  continuous  systems  this  is  defined  as: 

{Wei 0(f)  =  —BT  exp(—ATt)x 
For  discrete  systems  T1*  :  Cn  — >  /2[0,oo)  is  defined  as: 

(Kx)(t)  ee  BT(AT)k 
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Definition  4.4  For  a  given  a  state  space  representation  ( A,B,C ),  the  adjoint 
observability  operator  maps  past  outputs  to  internal  states  :  Cn  -*•  L2(~oo,  0] 

.  For  continuous  systems  this  is  defined  as: 

%y  =  [  exp(AT  t  )CTy{r  )dr 

J  —  OO 

For  discrete  systems  lIf '  :  Cn  — ’  1,{ — co ■  0]  is  defined  as. 

($»(()  =  E(^)‘CT»(r) 

—  CO 

The  importance  of  these  definitions,  aside  from  proving  reachability  and  de¬ 
tectability,  is  in  the  construction  of  the  controllability  and  observability  gram- 
mians  and  the  Hankel  operator. 

Definition  4.5  Given  a  stable  system  (A,B,C,D)  (  i.e.  the  eigenvalues  of  A 
of  a  continuous  system  are  all  contained  in  the  left  hand  side  of  the  complex 
plane,  for  discrete  systems  all  the  eigenvalues  are  contained  within  the  interior 
of  the  unit  disc),  the  continuous  controllability  grammian  is  defined  as: 

P=  r  exp{At)BBTexp{AT)dt  =  VeV*e  (4.2.3) 

Jo 

and  the  discrete  controllability  grammian  is  defined  as. 

P  =  f)  AkBBT{AT)k  =  (4-2.4) 

k=0 

Likewise,  the  observability  grammian  for  continuous  systems  is  defined  as: 

0  =  T  cxp(ATt)CTCexp(A)dt  =  (4-2.5) 

Jo 


34 


and  the  discrete  controllability  grammian  is  defined  as: 

CO 

Q  =  '£(AT)kCTCAk  =  V*oqo  (4-2.6) 

k= 0 

For  a  stable  continuous  system  these  grammians,  P  and  Q,  will  solve  the 
following  equations  [16]. 

AP  +  PAT  +  BBt  =  0  (4.2.7) 

AtQ  +  QA  +  CTC  =  0  (4.2.8) 

and  for  stable  discrete  systems,  these  grammians  will  solve. 

P  -  APAt  -  BBt  =  0  (4.2.9) 

Q  -  AtQA  -CtC  =  0  (4.2.10) 

Definition  4.6  For  a  given  system  with  impulse  response  h(t)  (  h(k)  for  dis¬ 
crete  systems),  the  Hankel  operator  is  defined  as  the  operator  that  maps  past 
inputs  to  future  outputs,  V  :  L2{- oo,0]  ->  L2[0,oo).  This  operator  for  continu¬ 
ous  systems  is  expressed  as: 

re o 

(ru)(t)  =  /  h{t  +  s)u(s)ds 

J  0 

and  for  discrete  systems  T  :  l2(—co,0]  — +  ?2[0,oo)  is  expressed  as. 

OO 

{Tu){n)  =  h(k  +  n)u(k) 
o 

For  linear  systems,  the  Hankel  operator  is  easily  constructed  using  the  con¬ 
trollability  and  observability  operators. 

(ru)(f)  =  Wc  (4.2.11) 
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4.2.2  Hankel  Operator  Properties 

The  fundamental  relation  between  the  Hankel  operator  and  the  controllability 
and  observability  grammians  is  shown  in  the  following  theorem. 

Theorem  4.1  The  spectral  decomposition  of  the  Hankel  operator  is  related  to 
the  eigenvalue  decomposition  of  PQ,  as  follows: 

{<rr}  =  {A'\PQ)} 

where  or  represent  the  singular  values  of  the  Hankel  operator  (  ref:  Francis  [9] 
(ch  5)). 

Proof:  let  X  be  a  nonzero  eigenvalue  o/PT  with  respective  eigenvector  v  then 

T“Tv  =  =  Xv  (4.2.12) 

Premultiplying  (  4.2.12)  by  %c  and  using  equations  (  4-2-3)  and  (  4.2.5)  and 
letting  a ;  =  gives: 

PQx  =  Xx  (4.2.13) 

Likewise,  premultiplying  equation  (4-2.13)  by  <STCQ  and  defining  u  =  **Qx  gives 
equation  (  4-2.12)  back  again. 

In  the  following  section,  the  performance  bounds  of  balanced  truncations  are 
shown.  These  performance  bounds  are  directly  related  to  the  singular  values  of 

the  Hankel  operator. 

All  identification  algorithms  for  this  structural  identification  experiment  are 
performed  within  the  discrete  domain.  The  desire  is  to  put  the  Markov  param- 
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eters  into  a  format  that  allows  for  easy  model  construction  and  takes  advantage 
of  the  truncation  properties  of  the  Hankel  singular  values.  This  is  demonstrated 
by  defining  the  controllability  and  observability  operators  in  a  matrix  form. 

Definition  4.7  Given  a  discrete  system  ( A,B,C )  and  input  u  6  l2(-co,0]  ex¬ 
pressed  as: 

u  =  [u(0)r,  it(— 1)T, . .  • ,  u(—n)T, . .  .]T 
The  infinite  controllability  matrix  operator  is  expressed  as: 

[  B  AB  A2B  ...  An~2B  An~1B  ...  ]  (4-2-14) 

Likewise ,  the  infinite  observability  matrix  operator  can  be  expressed  as: 

$o=  [  cT  ATCT  ( At)2Ct  ...  (AT)n~2CT  (AT)n~1CT  ...  ]T 

(4.2.15) 

Thus,  the  Hankel  operator  defined  by  equation  (  4.2.11)  is  in  fact  an  infinite 
dimensional  matrix  that  maps  h{-oc,  0]  ->  ^[0,oo)  with  the  same  singular 
values  as  shown  in  Theorem  4.1. 

The  advantage  of  the  discrete  Hankel  operator  is  that  it  can  be  constructed 
directly  from  the  Markov  parameters  or  impulse  response.  The  realization 
( A,B,C )  has  an  impulse  response  defined  by: 

h(k)  =  CAk~lB 
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The  infinite  dimension  matrix  T,  therefore,  can  easily  be  constructed  as  follows: 

h{l)  h(k  + 1)  ...  h(k  +  n)  ... 
r  =  h(k  + 1)  h{k  +  2)  ...  h(k  +  n  + 1)  ...  (4.2.16) 

C£  CAB  ...  CAfc_1B  ... 

CAB  CA2B  ...  CAfc  ...  (4.2.17) 

Equation  (  4.2.17)  is  identical  to  equation  (  4.2.11). 

Though  constructing  the  Hankel  matrix  is  simple,  several  factors  must  be 
considered. 

1.  The  Markov  parameters  h(k)  will  be  corrupted  by  noise.  For  the  experi¬ 
ment  of  this  report,  the  signal  to  noise  ratio  will  most  likely  be  low. 

2.  Computationally,  working  with  an  infinite  dimensional  matrix  is  impossi¬ 
ble.  Therefore,  the  Hankel  matrix  must  be  truncated  to  a  finite  order  of 
n.  The  effect  of  truncating  this  matrix  must  be  quantified. 

3.  The  experiment’s  goal  is  parameter  identification,  thus  the  effects  of  noise 
corruption  must  be  quantified. 

The  next  section  presents  an  error  bound  on  model  truncation  and  provides  an 
explanation  of  the  noise  effects  on  the  parameter  identification.  Following  the 
next  section,  the  effects  of  truncating  of  the  infinite  matrix  are  investigated. 
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4.2.3  Balanced  Realizations  and  Error  Bounds 

The  most  important  property  of  the  Hankel  operator  is  the  relation  between 
the  hankel  singular  values  and  the  frequency  response.  This  relation  will  be 
exhibited  with  the  error  bound  of  a  balanced  order  reduction. 

While  the  realization  (A,  B,  C)  is  invariant  to  transformations,  the  grammi- 
ans  P  and  Q  are  not.  Moore  [25]  was  the  first  to  appreciate  this  and  showed 
that  there  exists  a  transformation  such  that  the  grammians  will  become  diago¬ 
nal.  That  is: 

TPTt  =  diag{T,u  E2,0,0) 

{Tt)-1QT~1  =  diag{  EiAEa.O) 

where  the  matrices  E;  are  diagonal  matrices  with  elements  {c^,  a2, . . . ,  crni}  ar¬ 
ranged  in  descending  order  o\  >  ct2  >  • .  •  >  >  0-  The  transformed  represen¬ 

tation  (A,B,C)  =  ( TAT~\TB,CT -1)  will  become: 
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c=  Cc,o  0*0 

where  the  subscripts  denote  controllable,  uncontrollable,  observable,  and  unob¬ 
servable  subspaces. 

Given  this  transformation,  the  controllable  and  observable  portion  of  the 
representation  will  satisfy  the  following  Lyaponov  equations  (for  continuous  sys¬ 
tems). 

Ac,0Ei  +  £  !Aj0  +  Be,0Bl0  =  0 
+  EiAC)0  +  Cj0Cc,o  —  0 

Considering  only  the  minimal  portion  of  a  representation,  balanced  model  trun¬ 
cations  are  performed. 

Definition  4.8  Given  an  nth  order  minimal  balanced  representation  ( A ,  B,  C ,  D), 
with  P  =  Q  =  diag(Bi,  £2),  where  Sr  =  diag{erlt . . .  ,0*)  and  £2  =  dzag(<rk+1, . . . ,  crn) 
arranged  in  descending  order  o’!  >  cr2  >  . . .  >  crn  >  0,  and  the  realization 
( A,B,C )  conformally  partitioned  as: 
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Then  a  kth  order  balanced  reduced  representation  is  (An,  Hi,  Ci,  T)).  (Reference 
[24]  p.398) 

Using  the  definition  of  the  transfer  function  of  G(s)  =  D  +  C(sl-  A)_1B,  it 
was  proven  independently  by  Glover (19S4)  and  Enns  [8]  using  optimal  Hankel 
norm  approximations  that  Lx  norm  of  the  difference  between  the  nih  order  true 
system  Gn(s)  and  the  internally  balanced  truncated  kth  order  model  Gk(s)  has 
the  following  bound. 

n 

\\G{s)n  -  Gfc(s)||oo  =  sup  \Gn{jw)  -  Gk(jw)\  <  2  X)  (4.2.18) 

"  w  we(-oo,co)  i=k+ 1 

The  following  engineering  results  may  be  deduced  from  the  balanced  order 
model  reductions. 

1.  The  spirit  of  balanced  order  truncations  is  that  the  least  observable  and 
controllable  modes  are  the  ones  being  truncated.  The  closer  to  zero  the 
truncation  singular  values,  the  less  observable  or  controllable.  In  this  sense 
the  truncation  makes  intuitive  sense,  since  these  modes  are  of  least  concern 
to  a  controller  designer. 

2.  The  Too  bound  is  significant  if  ak  >  °k+ 1  and  ak+ 1  is  small.  This  implies 
that  the  sum  of  the  tail  of  singular  values  is  small  and  thus  the  Too  bound 

is  tight. 

3.  The  truncations  can  perform  well,  even  when  ak+ 1  is  close  to  ak.  Glover 
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(1984)  showed  this  could  be  true,  since  as  ak+i  — >  P°les  °f  trun¬ 

cation  can  approach  the  imaginary  axis.  In  fact,  as  a  pole  tends  to  the 
imaginary  axis,  it  tends  to  become  either  uncontrollable  or  unobservable. 
Thus,  the  truncated  modes  will  have  less  affect  on  the  model  s  perfor¬ 
mance.  (Theorem  4.2  [11]). 

4.  If  the  realization  (A,  B,  C )  was  constructed  from  noisy  Markov  parameters 
(ie.  ERA  realization),  and  if  an  estimate  of  the  noise  variance  exists,  then 
a  comparison  of  the  bound  and  the  noise  power  spectrum  is  useful. 
Given  a  noise  power  spectrum: 

oo 

$noise(w)  =  Rnoise{T)e 
—  cc 

if  the  Loo  bound  satisfies: 

n 

2  V)  (7i  ~  sup  I  (u>)  i 

»=fc+ 1  w 

then  the  effects  of  the  truncation  is  limited  to  the  removal  of  modes  that 
model  the  noise  contributions. 

4.2.4  Optimal  Hankel  Norm  Model  Reductions 

An  alternative  model  reduction  scheme  exists  that  provide  the  same  Loo  bounds. 
This  approach  is  the  optimal  Hankel  norm  model  reduction  and  was  first  pre¬ 
sented  by  Glover  (1984)  using  balanced  realizations.  The  advantage  of  this 
approach  is  that  in  addition  to  the  bound  the  Hankel-norm  bound  is  mini¬ 
mized.  Begin  first  with  the  definition  of  the  Hankel  norm. 
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Definition  4.9  The  Hankd-norm  of  a  system  G(s)  with  realization  (A,B,C) 
is  the  largest  singular  value  of  its  Hankel  operator  and  is  defined  as. 

||<S'(s)||//  =  *(rG)  =  ^max(PQ)  (4.2.19) 


Glover  showed  that  the  Hankel  the  norm  gives  the  L 2  gain  for  all  past  outputs 
to  all  future  outputs. 


\\G(s)\\h  =  sup 


IIj/||l2(o,oo) 


(4.2.20) 


x6L2(-oo,0)  IM|l2(-oo,0) 

The  advantage  of  the  optimal  Hankel  norm  model  reduction  is  that  the  Hankel 

norm  of  the  difference  between  the  true  model  and  the  reduced  k  order  estimate 

is  minimized  and  has  the  following  Hankel  norm. 


II G(s)  -  G(s)\\h  =  <?*:+ 1 

Computationally,  the  calculation  of  the  optimal  Hankel  norm  model  reduction 
is  performed  without  a  balancing  transformation  [22].  Nevertheless,  for  our 
experiment  the  balanced  order  reduction  approach  is  usually  chosen  since  it  can 
be  accomplished  during  construction  of  the  realization  (See  section  on  ERA). 


4.3  Modeling  Infinite  Dimensional  Systems  with  Finite  Dimensions 

Equation  (  4.2.17)  shows  that  the  discrete  time  Hankel  operator  can  be  con¬ 
structed  exactly  from  the  Markov  parameters  (  i.e.  the  impulse  response  sam¬ 
ples).  Once  constructed,  the  Hankel  operator  may  be  used  to  construct  a  re¬ 
alization  (A,B,C)  that  will  exactly  duplicate  the  system’s  impulse  response. 
Properties  of  this  linear  realization  are. 
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1.  An  nth  order  realization  will  exactly  match  the  first  2 n  Markov  parameters. 
The  model’s  Markov  parameters  later  than  2n  will  differ  from  the  true 
systems  impulse  response  due  to  noise,  nonlinear  components,  and  higher 

order  modes. 

2.  It  requires  an  infinite  dimensional  system  to  exactly  represent  the  true 
system’s  impulse  response  with  added  noise. 

The  second  constraint  can  be  visualized  by  looking  at  the  frequency  char¬ 
acteristic  of  white  noise  (  See  Figure  4.1  ).  In  order  to  exactly  represent  that 
frequency  response,  the  model  would  require  an  infinite  number  of  poles  inside 

the  unit  disc. 

Unfortunately,  it  is  computationally  impractical  to  consider  handling  infinite 
dimensional  matrices.  Additionally,  the  experimental  conditions  preclude  any¬ 
thing  but  a  finite  test  window.  However,  error  bounds  resulting  from  truncating 
infinite  systems  to  finite  systems  have  been  investigated  by  Curtain  et  al  [14]. 
The  error  bounds  determined  by  Curtain  will  now  be  presented. 

In  order  to  consider  infinite  dimension  space  the  following  definition  is  re¬ 
quired  [3] . 

Definition  4.10  An  infinite  dimensional  discrete-time  realization  of  a  system 
( A,B,C )  operating  on  a  real  separable  Hilbert  space  H  is  as  follows: 

—  Axn  4“  SlAn 

yn  =  Cxn  +  Dun 
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Figure  4.1:  Modelling  White  Noise  with  Poles  inside  the  Unit  Disc 
and  A,  B,  C,  D  are  linear  operators  such  that: 

A:  H  ^  H  B  :&m  H  C  :  H  -+W  D  :?Rm 

The  transfer  function  described  by  this  infinite  dimensions  model  G{s)  is  defined 
by  : 

G{z)  =  C{zl  -  A)~1B  +  D 

4.3.1  Measuring  the  distance  between  truncated  systems  and  infi¬ 
nite  systems 

In  order  to  measure  the  effects  of  truncating  an  infinite  Hankel  operator  the 
following  norms  will  be  used: 
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Definition  4.11  A  system  is  power-stable  if  there  exists  M  >  0  and  r  £  [0,1) 
such  that  ||Afcx||  <  Mrk  for  all  x  £  H  and  k  £  {0,1,2,...}  (Ref  [3]).  The 
Hankel  operator  T  associated  with  a  power-stable  system  is  nuclear  and  satisfies 

the  following: 

OO 

^  °k  <  co  (4.3.1) 

fc= i 

The  summation  in  equation  (  4-3-1)  is  referred  to  as  the  nuclear  norm,  || - 1| at- 
All  bounded  finite  systems  are  nuclear,  however,  this  is  not  strictly  true  with 
infinite  systems.  Note  that  this  norm  is  the  trace  of  the  singular  values,  thus  it 
is  also  called  the  trace  norm  [26]. 

Definition  4.12  A  power  bounded  or  l2-stable  system  will  satisfy  the  following 
constraint: 

(f>£)1/2<°°  (4-3-2) 

fc= i 

The  summation  in  equation  (4-3.2)  is  referred  to  as  the  Frobenius  norm,  ||.||f- 

The  nuclear,  Frobenius,  and  Hankel  norms  all  satisfy  the  basic  norm  prop¬ 
erties  of: 

1.  imi  =  q-||a|| 

2.  ||  A||  >  0  with  equality  if  and  only  if  A  =  0 

3.  [I A  +  B||  <  11  AH  +  ||jB||  The  triangle  inequality 
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Additionally,  these  norms  are  unitarily  invariant.  That  is  for  any  matrix  U  that 
satisfies  UUT  =  I  the  norm  of  T  is  unaffected  by  multiplication  with  U. 

||LT||  =  \\TU\\  =  ||r|| 

4.3.2  Bounds  and  Convergence  of  Truncated  Systems 

Systems  with  the  property  of  nuclearity  have  an  important  bound  on  the  h  norm 
of  the  impulse  response.  This  bound  was  shown  by  Curtain  (1988) [Theorem  2.1] 
and  the  bound  for  discrete  systems  is  stated  without  proof. 

OO 

s>(0l  =  w.  <2||r||«  (4-3-3) 

t=0 

What  must  be  quantified  are  the  effects  of  truncating  the  infinite  Hankel  operator 
matrix  in  equation  4.2.11  to  a  finite  order  matrix.  Let  Tn  represent  the  nth 
order  approximate  to  the  complete  Hankel  operator  T.  Of  interest  is  the  size  of 
appropriate  norms  on  the  difference  between  the  truncated  and  infinite  operators 

lir-rj. 

Curtain  et  al  showed  that  if  the  system  impulse  response  h(t)  €  Xi  0  T2  and 
the  Hankel  operator  was  nuclear,  then  the  truncation  of  the  Hankel  operator 
will  converge  to  the  true  Hankel  operator  as  the  order  n  gets  large.  That  is: 


||rn-r||w-o 

(4.3.4) 

II An  -  A||!  -  o 

(4.3.5) 

II  An  -  A || 2  -  0 

(4.3.6) 
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These  intuitive  conclusions  for  a  stable  system  also  come  with  error  bounds 
on  performance.  This  error  bound  is  similar  to  the  one  for  balanced  order 
truncations.  The  bound  is  as  follows: 

||G-Gn||oo  <2£>« 

i>n 

Where  G  and  Gn  are  the  transfer  functions  associated  with  the  respective  Hankel 
operators.  Additionally,  there  exist  performance  bounds  on  the  Lx  and  L2  norms 
on  the  difference  between  the  truncated  and  infinite  impulse  responses. 

4.3.3  LACE  nuclearity  assumptions  and  measurement  noise 

A  necessary  assumption  of  Curtain  et  al  (1988),  is  that  the  system  is  nuclear 
and  an  element  of  h  fl  L2  functions.  Nuclearity  is  dependent  on  the  system 
under  study.  It  was  shown  by  de  Vries  [7]  that  a  system  constructed  via  the 
Laplace  transform  of  a  partial  differential  equation  model  of  a  deployed  space¬ 
craft  essentially  equivalent  to  the  LACE  spacecraft  is  in  fact  nuclear.  A  neces¬ 
sary  assumption' is  that  the  booms  have  viscous  damping  (See  chapter  2).  This 
assumption  is  quite  reasonable  since  the  structural  damping  factor  tested  on 
ground  was  greater  than  1  percent.  Additionally,  any  damped  structural  system 
.  is  definitely  an  element  of  Lx.  Thus,  the  assumption  of  nuclearity  is  reasonable. 
Unfortunately,  noisy  experimental  data  will  not  be  nuclear.  Since  the  all 
measurements  are  corrupted  by  noise,  the  assumption  that  the  measured  system 
is  an  element  of  h  will  not  be  valid.  Referring  to  equation  (  4.3.3),  it  is  clear 
that  the  assumption  of  nuclearity  of  measured  data  will  be  violated. 
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Thus,  any  model  that  includes  the  effects  of  measured  noise  will  nullify  pos¬ 
sible  bounds  on  the  truncation  of  the  Hankel  operator.  While  this  might  seem 
disappointing,  it  certainly  is  a  logical  outcome.  However,  modeling  that  is  done 
in  a  stochastic  sense  will  not  be  confined  to  attempting  to  model  noise.  The  con¬ 
vergence  of  equations  (  4.3.4)  through  (  4.3.6)  may  indeed  occur.  This  difference 
between  deterministic  modeling  with  truncation  and  stochastic  modeling  with 
truncation  is  the  primary  advantage  of  the  approaches  MME  and  least-squares 
modeling.  ; 


4.3.4  Convergence  Criteria 


In  order  to  calculate  the  Nuclear,  Frobenius  and  Hankel  norms,  the  finite  op¬ 
erator  Tn  will  be  augmented  with  zeros  such  that  its  dimensions  match  that  of 
the  complete  operator. 


rn  = 


rn  o 

0  0 


norm 


The  trace  norm  has  the  property  that  for  balanced  hankel  operators,  the 
L-m  of  ||r  -  r„||jv  will  equal  the  difference  of  the  independent  trace  norms. 


||r-rn||/v  =  urn*  -  ||r„||iv 


The  Hankel  and  Frobenius  norms  do  not  have  this  property,  therefore,  the  tri¬ 
angle  inequality  must  be  used  to  bound  the  performance  of  these  two  norms. 
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How  a  truncated  Hankel  operator  characteristics  converges  to  the  character¬ 
istics  of  the  infinite  Hankel  operator  will  be  measured  using  a  Cauchy  sequence 
of  operators  T,  where  i  :  1  -»  oo.  The  quality  of  convergence  will  be  determined 

how  close  Ti  is  to  F1+i  that  is: 

||r,+i  -  Till  <  e 

if  e  is  sufficiently  small  and  is  converging  toward  zero,  then  the  confidence  in 
rn  with  n  >  i  will  be  high.  The  trace  and  Hankel  norms  will  be  of  particular 
interest.  Note:  The  Hankel  norm  represents  the  gain  of  the  system  (  Hi  gain  ). 

4.4  Eigensystem  Realization  compared  to  Balanced  Order  Reduc¬ 
tion 

The  most  remarkable  improvement  in  identification  schemes  was  the  Eigensys¬ 
tem  Realization  Algorithm  (ERA)  proposed  by  Juang  and  Pappa  (1985).  The 
idea  behind  the  algorithm  is  the  realization  algorithm  of  Kalman-Ho  [2],  but 
it  utilizes  the  numerically  robust  approach  of  singular  value  decomposition.  In 
this  section  it  will  be  shown  that  ERA  with  its  model  reduction  based  on  sin¬ 
gular  values  is  functionally  equivalent  to  the  Schur  method  of  balanced  order 

reduction  presented  in  1989  by  Safonov  et  al  [21]. 

The  Schur  method  is  based  on  the  large  singular  values  of  the  Hankel  operator 
and  its  singular  value  decomposition.  Singular  value  decomposition  is  defined 
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as  follows  (See  Maciejowski  [24]): 


H  =  U  ZVT 


(4.4.1) 


U  =  is  the  matrix  of  the  left  eigenvectors  of  HHT .  That  is 
HHT  =  U'E2UT. 

VT  =  is  the  matrix  of  the  right  eigenvectors  of  HT H.  That  is 
HtH  =  VZ2VT. 

E  =  diag[ai,  cr?, . . . .  <7;.]  where  k  equals  the  rank  of  H. 
ff,-  =  A  1'2{HHT) 

Safonov  et  al  (1989)  used  singular  value  decomposition  to  develop  an  algo¬ 
rithm  to  perform  order  reduction  without  requiring  balancing  the  system  first. 
(Note:  balancing  transformations  such  as  proposed  by  Moore  (1981)  can  be 
extremely  numerically  unstable).  The  approach  performed  an  eigenvalue  de¬ 
composition  of  the  matrix  PQ.  The  matrices  Vl,big  and  Vr ,big  composed  of 
the  right  and  left  eigenvectors  associated  with  the  “big”  eigenvalues  are  used  to 
perform  model  reduction. 
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4.4.1  Eigensystem  Realization  Algorithm  Construction  compared 
to  Schur  method  Order  Reduction 

The  ERA  approach  assumes  that  the  system  response  sampled  is  produced  by 
a  linear  system’s  impulse  response.  First  assume  a  single  input  multi-output 
(SIMO)  system  (A,B,C).  Working  with  this  realization,  the  discrete  system 
impulse  response  is  described  by: 

Y(k)  =  CAk~1B 


If  the  response  is  considered  a  free-decay  from  an  initial  condition  x(0)  the 
response  will  be  as  follows: 


Y(k)  =  CAk~lx(  0) 


From  here  forward,  the  free-decay  case  initial  condition  x(0)  will  be  equated 
with  input  matrix  B. 

Now  consider  H(k)  to  be  the  Hankel  operator  time  shifted  by  k  time  units. 


Y(k)  Y(k  +  1)  ...  Y(k  +  n) 


H(k )  = 


Y(k+  1)  Y{k  +  2) 


Y(k  +  n  +  l)  . . . 


It  is  easy  to  show  that  this  matrix  can  be  constructed  by  using  the  controllability 
and  observability  operators,  Tc  and  (see  equations  (  4.2.11)  and  (  4.2.17)). 


H(k)  = 
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It  can  be  shown  that  the  zero  delayed  Hankel  operator  H( 0)  has  a  left-pseudo 
inverse: 

H{  0)#  =  (i/(0)TJff(0))-1^(0)  (4.4.2) 

which  satisfies  the  following  condition: 

VcH(0fqo  =  I 

I  is  the  identity  matrix  of  appropriate  size.  Now  using  the  same  notation  as 
equation  (  4.4.1),  the  zero  delayed  Hankel  operator  has  a  singular  value  decom¬ 
position  of: 

H(Q)  =  UEVT  (4.4.3) 

'  H(  0)#  =  VZ-1UT  (4.4.4) 

Using  equation  (  4.4.2),  the  realization  (A,  B,  C )  is  constructed  from  the  Hankel 


operator  as  follows  (See  Ref  [13]). 

Y(k  + 1)  =  EpCAkBEm  (4.4.5) 

=  -  EjH{k)Em  =  E^0AkycEm  (4.4.6) 

=  EjqoycH*{0)yoAkycH{0)*yoycEm  (4.4.7) 

=  EZH{0)V(2-1UTVoAkVeVZ-1UTH(Q)En  (4.4.S) 

=  E%H{0)V('2-1UTH{l)V)kJr1UH(0)Em  (4.4.9) 


=  E*UZ{1/2\E(~1/2)UT H{l)Vrt~1/2))kZ{l/2)VT Em  (4.4.10) 

Ej  =  [7p, 0,0,  •  •  ■  ,]r  and  E £  =  [1, 0, 0,  •  •  •  ,]r.  A  full  order  realization  is 
constructed  with  A  =  Tj~1Ut H(1)V ,  B  —  Em,  and  C  =  EjU^1^2^.  A 
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realization  that  reduces  order  in  a  balanced  sense  will  select  only  the  “large” 
singular  values  of  the  Hankel  operator  in  the  construction.  Letting  £/.  = 
diag(ai,  •  ■  •  ,crk)  representing  the  kth  largest  singular  values  and  using  matrices 
Uk  and  14  that  select  the  first  k  columns  of  U  and  V,  the  construction  becomes: 

A  =  Sl_1/%Jtf(l)14-2jf1/2) 

B  =  ^/2)V?Em 
C  =  E*UkE  i1/2) 

4.5  Improvements  to  Eigensystem  Realization  using  Minimum  Model 
Error  Estimation 

While  ERA/Balanced  order  reduction  and  optimal  Hankel-norm  model  reduc¬ 
tion  schemes  provide  nice  error  bounds  in  the  frequency  domain,  they  do  not 
provide  any  bounds  or  insight  to  what  is  happening  in  the  time  domain.  These 
realization  approaches  are  deterministic,  exactly  matching  2 n  impulse  response 
samples  with  an  nth  order  system.  An  approach  presented  by  Mook  et  al  [4], 

[5]  addresses  time  domain  considerations.  This  section  explains  this  technique 
and  applies  it  in  a  stochastic  sense.  The  technique  coupled  with  a  realization 
and  truncation  scheme  (  such  as  Balanced  order  or  Optimal  Hankel-norm  model 
reduction)  will  produce  results  that  model  the  system  in  the  time  domain  with 
a  desired  error  variance  and  provide  improved  error  bounds  in  the  frequency 
domain.  The  necessary  background  is  presented  first,  and  then  followed  by  the 
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actual  minimization  scheme. 


4.5.1  Minimization  of  the  Objective  function 
Given  a  differential  system  that  is  modeled  as: 


x(t)  = 


(4.5.1) 


with  initial  state  x(t0)  and  a  cost  function  $(•)  associated  with  its  final  state 
x(tf).  Assume  a  cost,  L(u(-),  x(-),  •),  is  incurred  for  using  a  state  driving  input 
(note:  The  cost  is  a  function  of  input,  state  and  time).  The  integrated  cost  of 
this  input  will  be  ///  L(u(t),  x(t),  t)dt.  A  cost  function  J  is  then  constructed  to 
sum  all  related  costs.  This  cost  function  is  dependent  on  the  exogenous  input 
u(t)  and  is  defined  as: 


min  J  = 

U 


d[z(*/)]  +  /  '  L[x(t),u(t),t](It 

Jto 


It  was  shown  by  Pontryagin  that  this  system  can  be  optimally  driven  from  an 
initial  state  to  a  final  state  [17].  Pontryagin’s  optimization  principle  showed 
that  by  constructing  the  co-state  p(t)  defined  by: 

m—fam-j p  («-2) 

with  the  boundary  conditions: 


P(tf)  = 


df 

dx 


x[tQ )  specified 
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the  optimal  input  can  be  calculated  to  drive  x(t0 )  — ►  x(tf)  such  that  J  is  mini¬ 
mized.  The  optimizing  input  u(-)  satisfies 

dL[x(t),u{t),t]T  df(x(t),u(t),t)T }  3 

du  du 

Using  Pontryagin’s  optimization  principle  as  our  motivation,  a  system  with 
output  y(t)  is  used  to  optimally  model  the  true  outputs  y(t).  An  alternative 
cost  function  must  be  constructed  to  perform  an  optimization  that  accounts  for 
the  model  errors  at  each  time  sample.  The  cost  function  to  be  minimized  is 
constructed  as  follows: 


M  rtf 

min  J  =  T(yk  -  VkfR^iVk  ~Vk)  +  uT(r)lUu(r)dr  (4.5.4) 

u  k=i  Jt° 


yk  =  tk  sample  of  the  true  free-decay 

yk  =  tk  output  of  the  modeling  system,  where 

y(-)  =  $(*(•)>  «(*)»(•)) 

Rk  =  Noise  covariance  matrix  at  time  instant  k 

u(t)  =  input  to  system,  the  independent  variable  of  optimization 

W  =  Weighting  function  applied  to  inputs 

This  minimization  of  this  cost  function  has  two  goals.  The  first  is  to  minimize 
the  weighted  squared  difference  between  the  model  estimated  output,  the  true 
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system  output,  and  each  sample  instant.  The  weighting  function  Rk  is  the  noise 
covariance  of  the  true  system.  Normally,'  this  covariance  is  considered  to  be 
constant,  however,  it  may  vary  with  time.  Varying  the  noise  covariance  matrix 
is  very  useful  when  the  magnitude  of  the  noise  varies.  (The  noise  covariance 
variation  feature  will  be  utilized  in  the  analysis).  The  second  goal  is  to  minimize 
the  weighted  square  of  the  exogenous  input  u(t).  This  weighting  function  may 
be  varied  as  well.  While  there  is  no  advantage  in  varying  its  value  over  time,  its 
magnitude  is  important  in  our  experiment. 

These  two  goals  are  competing  within  the  minimization  process.  The  opti¬ 
mization  result  depends  on  the  importance  given  the  respective  weighting  fac¬ 
tors.  For  example,  if  the  input  weighting  factor  W  is  kept  small,  then  the 
optimization  will  calculate  an  input  that  will  cause  the  system  model’s  output 
to  near  exactly  track  the  true  output  at  each  sample  instant  tk.  If  the  input 
weighting  factor  is  large  ,  then  the  optimization  will  calculate  an  optimal  input 
that  contributes  very  little  change  to  the  modeled  system’s  output. 

Given  this  background,  an  improved  approach  is  now  outlined  (Ref  [5]): 

1.  Optimally  calculate  an  input  u(t)  so  that  the  error  covariance  matrix  is 
approximately  equal  to  the  true  noise  covariance  matrix  as  follows  (Ref 

[4]): 

Rll „(0)  =  -1  E(w  -  yt)T(ytc  -  h)  ~  (4-5-5) 

M  Jfc=l 

2.  Calculate  the  modeled  output  yk  using  the  optimizing  input  u(i). 
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3.  Calculate  a  realization  using  the  modeling  output  yk  instead  of  the  true 


measurements  yk- 


4.  Calculate  the  error  bounds  of  this  realization. 


4.5.2  Optimization  Problem  formulation 


The  calculation  of  the  optimizing  input  u(-)  will  require  solving  differential  equa¬ 
tions  with  boundary  values  for  the  states  and  the  co-states.  Equation  (  4.5.4) 
will  be  re  written  as:  - ' 


min  J  = 

U 


M  rtf 

'£Kl(x(ti),ti)+  uT(r)Wu(r)d 

i= l  Jt° 


(4.5.6) 


The  optimized  co-states  must  include  jumps  discontinuities  to  account  for 
minimizing  the  cost  of  the  error  at  each  measurement  (Ref  [4]  and  [12]).  This 
discontinuity  results  from  the  summation  term  in  equation  (  4.5.6)  and  requires 
that  the  co-states  be  adjusted  at  each  output  sample.  Thus,  if  at  time  tj  a 
sample  is  taken,  the  co-states  are  adjusted  by  the  derivative  of  the  summation. 

A((t)  =  A(fJ)  -  (4-5'7> 


For  the  experiment  under  consideration,  the  modeling  system  will  be  chosen  to 
be  linear  and  described  by  equations  (  4.2.1)  and  (  4.2.2).  Using  this  assumption, 
equations  (  4.5.1)  through  (  4.5.7)  take  on  a  much  simpler  form. 


A  =  -At  A 

(4.5.8) 

2  WTu{t)  =  -Bt  A 

(4.5.9) 
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A  (it)  =  \{tj)  -  2  CTRkl(y{tk)  ~  Cx{tk))  (4.5.10) 

Note:  It  can  be  shown  that  all  optimizing  inputs  u(-)  are  in  range  of  B  . 
Mook  [6]  suggested  a  problem  formulation  of  the  problem  so  that  it  would 
lend  itself  to  a  linear  algebra  solution.  This  approach  takes  the  n  internal  states 
and  the  n  co-states  and  combines  the  two  parts  into  a  2 n  order  system.  This 
new  system  is  a  first  order  differential  problem  with  initial  and  final  boundary 
conditions.  The  2 n  system  will  be  described  as. 

z{t)  =  Az  (4.5.11) 

where:  z(t)  =  [x\(t),X2(t), . . .  ,xn(t),pi(t),p2(t), .  ■ .  ,Pn(t)] 

The  2 n  x  2n  matrix  A  is  constructed  from  the  original  linear  differential 

equations. 

A  -\B(WT)~lBT 
A  = 

0  -AT 

The  differential  equation  (  4.5.11)  satisfies  boundary  values  at  t0  and  tj.  The 
boundary  conditions  for  the  co-states  are  assumed  to  be  zero  at  the  beginning 
and  the  end,  p(t0)  =  p{tf)  =  0.  This  assumption  is  reasonable  since  there  is 
no  information  about  the  internal  states  x.  [Note:  If  the  free-decay  assumed  to 
be  due  to  an  initial  condition,  then  the  initial  of  x(tc)  =  B  may  be  assumed.] 
Using  the  assumption  on  co-state  boundaries,  the  boundary  value  equation  is 

constructed: 

0  =  B0z(t0)  +  BfZ(tf)  (4-5.13) 


(4.5.12) 


59 


where  0  =  \pi(t0),p2(t0), . . .  ,pB(io),Pi(*/), !*(*/), .  -  .  ,pn(tf)]T  =  [0r,0T]r  and 


Ba  = 


Bf  = 


0  I 
0  0 

0  0 
0  I 


where  all  submatrices  are  dimension  n  x  n. 

Finally,  the  jump  discontinuities  for  the  2 n  system  are  expressed  as: 

r(t+)  =  Dthz(tj;)  +  Yth 


(4.5.14) 


where  yk  is  the  kth  sample  of  the  true  system  and  Dtk  and  Ytk  are  defined  as: 


I  0 

-2  CTR;lC  I 


Yt,.  = 


0 


2Ctr;1 


4.5.3  Solving  the  Two-point  boundary  value  problem  by  Multiple 


Shooting 


One  approach  to  solving  the  two-point  boundary  value  problem  (TPBVP)  is  the 
approach  of  shooting,  That  is,  guess  at  the  form  of  the  solution  and  formulate 
the  problem  so  as  to  solve  the  solution. 


y(t)  =  v(t)  +  V(t)c 
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such  that  the  boundary  condition  equation  (  4.5.13)  is  satisfied: 

B0[v(t0)  +  V(t0)c\  +  Bf[v(tf)  +  V(tf)c }  =  (3  (4.5.15) 

where: 


1.  c  is  the  vector  of  coefficients  for  the  time  interval. 

2.  v(t)  is  a  particular  solution  2 n  vector.  The  particular  solution  will  be 
assumed  to  have  zero  initial  condition  v(t0)  =  0. 

3.  V{t)  is  a  solution  2 n  x  2n  matrix.  The  solution  matrix  will  be  assumed 
to  have  zero  initial  condition  V(t0 )  =  I. 

Thus  using  the  state  space  representation  of  equation  (  4.5.11),  the  solution  has 
a  simplified  form: 

y(tf)  =  exp(A(tf  -  t0))c 

Letting  Atf~t0  =  exp{A{tf-  ta))  and  using  the  initial  conditions  for  v{t)  and 
V(t )  the  solution  of  equation  (  4.5.15)  can  be  expressed  as. 

c  =  [B0  +  BfAtj-t.]-1?  (4-5-16) 

The  solution  to  c  will  exist  if  the  inverse  to  [Ba  +  BfAtf-t J-1  exists.  Mook 
[6]  showed  that  this  is  a  necessary  and  sufficient  condition  to  solve  the  simple 
shooting  problem.  However,  if  the  differential  equation  is  stiff,  then  the  solution 
may  end  up  numerically  unstable  [23]. 
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The  problem  of  stiff  differential  equations  can  be  solved  by  breaking  the 
time  interval  into  subintervals  where  the  boundary  conditions  between  adjacent 
subintervals  are  consistent.  If  m  subintervals  are  created,  then  there  will  be 
m  +  1  simultaneous  equations  in  the  form  of  equation  (  4.5.16)  to  solve  in  order 
to  get  the  coefficients  for  each  subinterval.  This  technique  is  called  multiple 

shooting. 

Multiple  jump  discontinuities  are  easily  encorporated  into  the  process.  As¬ 
suming  measurements  taken  at  time  instants  ft  with  k  €  {1,  •  •  • ,  M},  then 
during  any  shooting  time  interval  equation  (  4.5.14)  becomes: 

+  Ytk+ 1 

where  tfc,4*+i  are  within  shooting  interval. 

Thus  with  this  formulation,  the  solution  to  the  TPBVP  becomes  one  of  linear 
algebra  to  solve  for  the  coefficients  of  each  subinterval. 

4.6  Model  Identification  using  Least-Squares  model  fitting 

Previous  sections  have  presented  approaches  for  realization  of  a  system  assuming 
that  noisy  samples  of  the  impulse  response  are  given.  In  this  section  an  alternate 
identification  scheme  is  proposed  given  the  same  assumptions.  The  difference 
between  this  algorithm  and  the  previous  algorithms  is  that  the  identification 
will  be  done  entirely  in  a  stochastic  sense.  This  identification  scheme  constructs 
a  realization  that  minimizes  the  norm  between  the  modeled  system’s  impulse 
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response  and  the  true  measurements.  The  advantage  of  a  stochastic  based  real¬ 
ization  will  occur  when  the  number  of  measurements  becomes  relatively  high. 

The  realization  schemes  such  as  ERA  are  deterministic  and  exactly  dupli¬ 
cate  the  first  2 n  samples  of  a  system  with  a  nth  order  system.  While  this  may 
be  acceptable  for  sma-ll  values  of  n,  computationally,  it  becomes  weak  as  n  be¬ 
comes  large.  Additionally,  the  desired  error  bounds  become  less  significant  as 
the  number  of  modes  modeled  increases  due  to  the  effects  of  noise  (See  Sec¬ 
tion  on  Infinite  Dimensional  Systems).  This  stochastic  approach  will  use  the 
least-squares  approach.  [  Note:  Attempts  to  use  similar  techniques,  such  as 
instrumental  variables,  did  not  provide  any  improvement  and  were  more  likely 
to  produce  erroneous  overdetermined  systems  [1]]. 

4.6.1  Autoregressive  Models 

A  linear  dynamic  discrete  system  can  be  modeled  as  process. 

A(q)y{n)  —  B(q)u(n)  +  C(q)e(n) 

where: 

1.  q  is  the  delay  operator,  identical  to  z-1  in  the  z-transform. 

2.  u(n)  is  the  system  input  and  e(n)  is  unmeasurable  noise  input  to  the 
system. 

3.  A{q)  =  {I  +  Axq -1  +  A2q~2  +  •  •  •  +  Anaq~na ] 
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4.  B{q)  =  [I  +  B1fj~1  +  B2q~2  +  •  •  •  +  Bnbq-nb] 

5.  C(q)  —  [I  +  Ciq~l  +  C2q  2  +  •  •  •  +  Cncq  c] 

6.  (na,  rtf,,  nc)  are  the  orders  of  the  dynamics  of  the  different  parts  of  the 
model. 

The  transfer  function  of  this  system  written  in  matrix  fraction  description  is: 


y(n)  =  A(q)  1  B(q)u(n)  +  A(q)  1C(q)e(n) 

For  the  LACE  experiment,  white  noise  is  assumed.  This  causes  C(q)  to  be 
unity,  and  the  system  reduces  to  the  standard  ARX  (autoregressive-exogenous 
input)  format.  Also  assuming  the  system  is  in  free  decay,  u(n)  =  0,  the  system 
becomes  an  autoregressive  model.  The  one  step  ahead  predictor  of  y(n)  is  easily 
calculated  to  be  (Ref  [20]): 

y(n\9)  =  ipT{n)e 

where  0  is  the  matrix  of  coefficients  in  A{q )  and  <p(k)  is  a  column  vector  of 
system  samples: 

<?(k)  =  [-y{k),  ~y(k  -  1),  •  •  • ,  -y{k  -  na)]T 
9  =  [Ai,  A2,  •  ■  • ,  Ana]T 

The  error  associated  with  each  estimate  y(n)  is  then  the  innovation  e(n)  (Ref 

[20]): 

e(n)  =  y(n)  -  {n)9 
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4.6.2  Autoregressive  model  performance 


Letting  ZN  =  {y(n)|n  €  {1,  — ,  JV}},  the  autoregressive  approach  will  injec¬ 
tively  map  the  data  set  ZN  to  the  space  of  linear  model  of  desired  dimension. 

ZN  9 

The  goal  of  the  parameter  estimation  will  be  to  minimize  the  value  of  a  given 
cost.  The  cost  function  VN  associated  with  a  system  sampled  N  times  is  defined 
to  be:  • 

vN(e,zN)  =  ^jri(e(k,e)) 

n  k= i 

Using  the  least  squares  criteria  the  function  /(■)  is: 

Z(e(n))  =  ^er(n)A"1e(n) 

Where  A  is  the  weighting  function  on  the  error  components.  The  optimizing  set 
of  coefficients  is  calculated  by  linear  algebra  as: 

fc= 1  iV  tel 

Note:  That  by  letting  A  =  /,  the  first  matrix  in  equation  (  4.6.1)  becomes 
i  ¥>(fc)<Pr(fc)]  =  Rn(°)  the  N  samPle  correlation  matrix. 

The  LACE  experiment  the  cost  function  Vn{9,  Zn)  is  modified  to  account  for 
knowledge  of  the  noise  level.  This  is  accomplished  by  using  a  weighting  function 
fik.  The  cost  function  now  becomes: 

vk(0,zn)  = 

iV  k= 1 
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and  the  optimizing  least  squares  estimate  becomes: 

6„  =  (4.6-2) 

iV  k=  1  iV  *=1 

This  identification  algorithm  is  an  optimization  over  a  set  of  models  of  order 
na.  The  model  &is  is  then  an  element  of  the  set  of  linear  na  models  (  thetais  G 
Dna  C  Rna  for  SISO  systems  ).  This  optimization  is  guaranteed  to  converge  to 
the  true  model  90  as  N  ->  oo  provided  two  assumptions  hold. 

1.  e(n)  is  independent  white  noise. 

2  0  a  n  That  is,  the  true  model  is  indeed  an  element  of  the  set  of  models 
that  the  optimization  is  done  over. 

Assumption  (1)  is  reasonable,  given  the  nature  of  the  experiments  samples  (See 
chapter  1).  Assumption  (2)  however  is  not  guaranteed.  The  dynamics  under 
interest  are  generated  from  a  sampled  frequency  limited  infinite  dimension  sys¬ 
tem.  Thus,  while  a  finite  order  model  may  model  the  true  system  relatively 
closely,  it  cannot  be  assumed  that  90  G  DUa  when  na  is  finite.  Additionally, 
considering  the  potential  nonlinear  contributions,  the  case  form  assumption  (2) 
becomes  even  weaker. 

The  approach  chosen  for  this  experiment  is  to  use  a  sufficiently  large  value  of 
na  to  model  the  sampled  dynamics  to  within  a  desired  error.  Ljung  [19]  showed 
that  9is  will  converge  to  the  best  possible  solution  in  the  model  set  DUa  as  the 
number  of  samples  gets  large.  This  is  equivalent  to  minimizing  the  Z2[0,  N]  norm 
of  the  errors.  As  the  order  gets  large  the  tivo  things  occur. 
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1.  The  matrix  [jj  ZLi  Pk<p(k)<?T  {k)]  will  become  poorly  conditioned. 

2.  The  highest  orders  of  the  system  will  primarily  model  the  noise  and  non¬ 


linear  behavior. 

Since  the  goal  is  parametric  identification,  the  least  detectable  higher  or¬ 
der  modes  must  be  isolated  from  the  strongly  detectable  modes.  This  is  ac¬ 
complished  by  converting  the  autoregressive  model  into  a  impulse  response 
model  and  utilizing  the  balanced  order  or  optimal  Hankel  norm  order  truncation 
schemes.  The  procedure  is  as  follows: 

Step  1:  Calculate  d\s  of  an  order  such  that  Vn(0u,Zn )  is  sufficiently  small. 

Step  2:  Calculate  transfer  function  G(q)  =  A(q)~1B(q)  where  B(q)  is  chosen 
to  create  an  impulse  response  model. 

Step  3:  Truncate  system  to  dominant  modes  using  balanced  order  reduction  ( 
or  Optimal- Hankel  Norm  reduction). 

4.7  Sampling  Rate  consideration 

Computationally,  longer  data  sets  will  provide  better  identification  results.  How¬ 
ever,  oversampling  data  does  not  provide  any  advantages.  This  becomes  clear 
when  analyzing  the  resulting  parameter  identification  over  the  unit  circle.  An 
algorithm  can  only  guarantee  parametric  identification  performance  up  to  a  dis¬ 
tance  e  from  the  true  parameters.  As  the  data  ensemble  becomes  large,  a  good 
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identification  scheme’s  parameter  results  will  converge  to  the  true  values,  or  at 
least  sufficiently  close.  This  is  expressed  as: 

lim  e  — >  0 

•  N—*oo 

Each  parameter  then  is  guaranteed  to  be  close  to  its  true  value.  Unfortunately, 
as  modes  are  oversampled,  the  pole  placement  on  the  unit  disc  converges  near 
the  origin.  When  this  occurs,  there  are  two  effects: 

1.  Modal  frequency  identification  errors  f error  become  larger  for  a  given  e 

and  sampling  rate  cos.  That  is  f error  ~  While  us  — ►  Large,  the 

identification  convergence  (  es  — +  smaller)  will  not  necessarily  keep  up. 

2.  Multiple  modes  are  all  converging  to  the  origin,  causing  the  ball  Bc  of 
radius  e  associated  with  each  mode  to  overlap.  (  See  Figure  4.2). 

4  It  can  then  be  concluded  that  a  judicious  sample  rate  selection  will  separate 
the  most  observable  modes  from  each  other  and  decrease  the  error  in  estimated 
frequency.  If  many  modes  are  present,  filtering  and  frequency  selection  can  be 
used  to  best  isolate  and  identify  the  modes.  In  conclusion,  the  length  of  the 
data  ensemble  is  limited  by  period  of  observation  and  by  sample  rate  selection. 
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Figure  4.2:  Convergence  of  modes  on  unit  disc  due  to  oversampling 


4.8  Example  of  Identification  Technique  performance  using  LACE 
PDE  model  response 

An  example  of  the  performance  of  the  various  techniques  is  now  demonstrated, 
in  chapter  4  performance  on  actual  data  is  discussed  with  all  supporting  rational. 

The  data  used  was  generated  from  an  impulse  response  of  the  linear  PDE 
model  of  chapter  2.  An  internal  damping  factor  of  .01  was  used  (Note:  This  is 
less  than  1  percent  damping  at  frequencies  below  5  Hz.).  An  assumed  observa¬ 
tion  window  of  approximately  3  minutes  chosen  with  a  sampling  rate  of  5  Hz. 
The  peak-to-peak  excitation  was  set  to  approximately  10  mm/sec.  White  noise 
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with  a  covariance  of  1  mm/sec  was  added  to  the  impulse  response.  Additionally, 
bursts  of  noise  with  covariances  of  2  mm/sec  and  3  mm/sec  respectively  were 
added.  The  noise  covariance  information  is  assumed  to  be  known  a  priori  in 
the  identification.  Using  this  data,  identification  schemes  were  produced.  The 
true  modal  frequencies  were  solved  with  a  =  0  as  in  chapter  2.  In  Figure  4.3 
the  Fourier  transform  shows  that  not  all  modes  are  detectable  through  the  noise 
floor.  This  is  due  to  the  high  noise  level  and  the  relative  displacement  of  the 
boom  tip  for  the  modal  frequencies  (See  chapter  2). 

In  Table  4.1  the  results  of  the  various  identification  schemes  are  compared. 
The  ERA  approaches  used  only  50  samples  in  the  construction,  more  conver¬ 
gence  could  be  expected  by  using  more.  Autoregressive  technique  used  all  the 

data. 


true 

freq 

ERA 

ERA/MME 

AR  model 
truncation 

0.0191 

- 

- 

- 

0.1298 

0.1305 

0.1305 

0.1259 

0.2581 

- 

- 

- 

0.3238 

0.3275 

0.3275 

0.3174 

0.7567 

0.7679 

0.7637 

0.7642 

0.8217 

- 

- 

- 

0.8639 

0.S702 

0.8669 

0.8637 

0.942S 

- 

- 

- 

1.7061 

- 

- 

• 

1.7223 

- 

- 

* 

Table  4.1:  Modal  Identification  of  noisy  PDE  impulse  response 
In  Table  4.2  the  L ^  error  bounds  for  the  ERA  approach  and  the  stochastic 
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ERA/MME  approach  are  shown.  In  Figures  4.4  and  4.5  the  convergence  of 
the  singular  values  and  the  frequency  responses  of  the  approaches  are  shown. 
It  is  clear  that  the  MME  step  will  generate  a  nuclear  hankel  operator.  The 
covariance  chosen  in  the  MME  step  was  2.3,  which  is  approximately  equal  to 
the  true  average  noise  covariance. 


ERA 

ERA/MME 

440.5 

40.3 

Table  4.2:  bound  on  model  truncation 

Note:  This  example  exhibits  the  dependence  on  sampling  rate.  All  tech¬ 
niques  performed  adequately  for  identifying  the  higher  modes.  However  the 
lowest  modes  identification  could  be  improved.  A  second  iteration  done  with  a 
decreased  sampling  rate  will  provide  better  identification  of  these  modes. 
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Figure  4.4:  Singular  values:  ERA 
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Figure  4.7:  Frequency  Response:  ERA/MME 


Figure  4.8:  Frequency  Response:  AR 
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CHAPTER 

Five 

Analysis  of  Experiment  with  LACE  spacecraft 


In  early  January  1990  the  laser  Doppler  measurement  experiment  was  per¬ 
formed.  The  test  data  was  collected  at  Lincoln  Laboratories  and  the  data  was 
then  transferred  back  to  the  Naval  Research  Laboratories  for  analysis.  The  test 
coordinated  the  deployment  and  retraction  of  the  lead  boom  with  the  observable 
window  for  laser  illumination.  Four  observations  were  measured  during  the  test 
period.  Of  these  tests,  only  two  windows  observed  the  Doppler  returns  while 
the  spacecraft  experienced  free-decay  vibration.  Currently,  these  two  sample 
periods  are  the  sole  periods  used  in  the  structural  identification  process. 

The  first  observation  period  is  shown  in  Fig  5.1.  This  observation  win¬ 
dow  captures  the  dynamic  motion  of  the  lead  boom  while  the  boom  is  being 
retracted  and  immediately  following  the  end  of  retraction.  The  second  half  of 
this  observation  will  be  used  in  analysis.  The  complete  period  is  informative 
in  that  is  shows  that  the  dynamic  events  can  be  detected  reasonably  well.  At 
approximately  90  seconds  into  the  observation,  the  dynamic  motion  due  to  the 
retraction  of  the  boom  ends.  From  that  moment  on,  the  detected  motion  of 
the  boom  is  due  to  only  the  rigid  body  and  vibration  modes.  Later  the  rigid 


75 


body  motion  is  removed  by  calculating  the  predicted  motion  and  subtracting 

this  motion  from  the  observed  motion. 

The  second  observation  period  is  shown  in  Fig  5.2.  This  observation  window 
was  taken  approximately  10  seconds  after  the  retraction  was  completed.  Since 
there  is  some  delay  in  the  measurement  collection,  the  higher  modes  will  have 
decayed  more  in  this  window.  Fortunately,  this  observation  window  is  much 
longer.  The  increased  number  of  samples  makes  all  the  calculations  more  robust 
to  the  presence  of  noise.  Additionally,  the  longer  time  period  will  allow  a  more 
meaningful  analysis  can  be  done  when  decimating  in  time.  As  mentioned  in 
the  previous  chapter,  decimating  in  time  will  help  separate  the  lower  frequency 
modes  within  the  unit  disc. 

Hopefully,  in  the  future  the  opportunity  for  more  experiments  will  occur. 
As  shown  later  in  this  chapter,  the  evaluation  of  these  two  limited  observation 
windows  provide  very  promising  results. 
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5.1  Dynamic  Body  compensation 


In  addition  to  the  vibrational  motion,  there  is  the  rigid  body  motion  of  the 
spacecraft  is  due  to  the  change  in  the  ground  observer’s  aspect  angle.  This 
angular  change  is  shown  in  Figure  5.3.  The  apparent  change  has  two  effects  on 
the  observed  Doppler  shifted  laser  return. 

1.  The  rigid  body  motion  will  be  detected.  Referring  to  Figure  2.6  of  actual 
measured  data,  this  mode  roughly  approximates  a  mode  with  a  period 
100  to  200  seconds  or  .01  Hz  to  .005  Hz.  Thus  the  frequency  separation 
between  this  apparent  rigid  body  motion  and  the  lowest  vibrational  mode 
of  .019  Hz  is  very  small.  In  order  to  detect  the  lowest  mode,  this  apparent 
motion  must  be  eliminated. 

2.  The  observed  damping  factors  of  the  detected  modes  will  be  affected  by 
the  position  of  the  spacecraft  in  the  sky.  If  the  observation  window  is 
before  the  spacecraft  comes  overhead,  the  calculated  damping  factors  will 
be  smaller.  If  the  modal  frequency  is  small,  the  calculated  damping  factor 
could  indeed  become  negative.  Likewise,  if  the  observation  window  is  after 
the  spacecraft  goes  overhead,  the  calculated  damping  factor  will  become 
larger  than  the  true  value.  Figure  5.3  demonstrates  this  effect  by  showing 
the  projected  pitch  velocity  vectors  onto  the  line-of-sight  vector.  In  the 
figure,  the  velocity  vector  is  always  in  the  pitch  direction  with  respect  to 
the  spacecraft  and  is  kept  at  a  constant  level.  Thus,  the  projection  of  the 
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velocity  vector  onto  the  LOS  vector  increases  before  passing  overhead  and 
decreases  after  passing  overhead. 

The  effects  of  this  apparent  motion  needs  to  be  adjusted  for.  This  is  ac¬ 
complished  by  calculating  orbital  position  with  respect  to  the  observation  site 
and  calculating  the  apparent  motion.  The  calculated  rigid  body  motion  was 
simply  subtracted  from  the  observed  motion.  Then  the  change  in  aspect  angle 
was  corrected  for  by  dividing  by  the  cosine  of  the  angle  between  the  LOS  vector 
and  the  pitch  plane.  - ' 

5.2  Identification  Analysis 

The  in-phase  and  quadrature  data  measured  at  Lincoln  Laboratory  was  analyzed 
in  order  to  calculate  the  relative  motion  of  the  spacecraft.  The  Fourier  transform 
of  the  data  was  filtered  such  that  the  relative  difference  between  the  return  of 
the  lead  retro-reflector  and  the  bus  retro-reflector  would  the  track  the  dynamic 
motion  of  the  spacecraft. 

The  test  data  has  Doppler  shift  measurements  collected  at  a  62  Hz  sample 
rate.  This  rate  is  significantly  higher  than  the  modes  of  interest.  To  filter  out 
some  of  the  noise  present,  groups  of  5  samples  are  taken  and  the  median  value 
of  the  population  is  selected.  This  technique  reduces  the  sample  rate  to  12.4  Hz, 
which  is  more  than  adequate  for  analysis.  Reducing  the  sampling  rate  further 
will  be  done  with  filtering  and  decimation  in  time.  Using  a  larger  population 
and  selecting  the  median  will  provide  more  noise  rejection  than  averaging.  This 
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is  due  to  the  small  population  used  in  the  analysis.  Averaging  populations  of  5 
or  10  samples,  is  susceptible  to  noise  components,  thus  the  median  is  the  most 
effective.  Another  technique  is  to  choose  the  most  likely  Doppler  shift.  However, 
the  population  size  is  still  too  small  for  this. 

5.2.1  Identification  Results 

Using  this  data,  an  initial  analysis  was  performed.  All  the  data  was  initially 
filtered  to  eliminate  high  noise  components.  The  cutoff  frequency  was  chosen  to 
be  .75  Hz.  This  number  allows  for  decimation  in  time  and  provides  margin  for  the 
higher  modes.  The  data  sets  were  then  decimated  in  time  into  distinct  data  sets. 
Each  of  these  data  sets  were  then  analyzed  and  statistically  compared.  (Note. 
The  minimum  model  error  results  generated  optimal  data  such  the  variance 
between  the  actual  and  simulated  data  was  approximately  2.4  mm/sec.  This 
corresponds  to  one  bin  in  the  discrete  Fourier  transform  of  the  Doppler  shifted 
laser  return.)  Table  5.1  lists  information  about  each  observation  windows. 
Tables  5.2  through  5.5  details  the  results  of  this  study.  Figures  5.4  and  5.5 
plot  the  detected  frequencies  and  damping  factors  for  all  the  data  sets  of  both 

observation  windows. 

The  results  of  this  analysis  are  quite  impressive.  Given  such  short  observation 
windows,  the  algorithms  detected  modes  which  correlate  well  with  predictions. 
Of  the  first  four  vibration  modes,  three  were  detected  on  day  91010.  The  lowest 
mode  is  only  detectable  if  the  rigid  body  motion  is  removed  from  the  data.  One 
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of  the  eight  data  sets  for  day  91010  did  not  detect  the  lowest  mode.  Thus, 
a  asterisk  is  used  in  Tables  5.4  and  5.5  associated  with  the  lowest  mode  to 
designate  that  only  7  of  the  data  sets  was  used  in  the  statistics. 

On  day  91008,  only  two  of  the  four  were  detected.  The  day  91008  observation 
window  was  only  25  seconds  long  and  this  window  was  too  short  to  collect  enough 
data  on  the  lowest  mode  (modal  period  is  ~  50  seconds). 

Note,  there  is  one  moderately  damped  mode  ~  .52  Hz  which  is  detected  on 
both  days  that  the  neither  the  PDE  or  FEM  model  predicted.  This  mode  is 
probably  not  due  to  the  boom  elements  because  the  damping  is  too  high.  Ad¬ 
ditionally,  this  vibration  does  not  correlate  with  any  yaw  or  roll  mode  predicted 
by  analysis.  While  the  explanation  for  this  mode  is  yet  unknown  it  is  believed 
to  be  due  to  either  unmodeled  behavior  within  the  boom  canisters  or  vibrations 
that  are  due  to  the  cantilevered  retro-reflector  plate.  Since  the  retro-reflector 
center  of  mass  is  offset  from  the  boom  axis.  A  local  vibration  may  exist  at  the 
boom  tip  due  solely  to  the  retro- reflector  plate.  If  this  is  the  case,  the  vibration 
would  be  due  to  a  local  elastic  behavior  in  the  boom.  Such  elastic  behavior 
would  tend  to  be  nonlinear  and  will  have  higher  damping. 

The  nonlinear  vibration  assumption  is  supported  by  the  size  of  the  damping 
factor  for  the  .52  Hz  mode.  Its  damping  factor  was  between  5  and  11  percent. 
This  level  is  significantly  higher  than  the  other  modes.  This  difference  may 
be  attributed  to  the  nonlinear  behavior  in  either  the  boom  canister  or  at  the 
cantilevered  retro-reflector  plate. 
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The  identified  modes  were  within  .005  to  .01  Hz  of  the  predicted  modes 
depending  on  the  identification  technique  compared.  While  this  is  impressive, 
the  identification  of  the  damping  factors  had  a  larger  variance.  Of  the  predicted 
modes  that  were  identified,  Tables  5.3  and  5.5  show  that  two  of  the  three 
modes  had  damping  of  approximately  one  to  two  percent.  This  agrees  well  with 
experiments  done  on  the  ground.  The  damping  factor  for  the  lowest  mode  of 
.019  Hz,  had  a  large  variance  associated  with  it.  This  variance  is  due  to  the 
fact  that  the  period  of  observation  was  very  short  in  comparison  to  the  period 
of  the  mode.  In  the  future,  longer  observation  windows  will  provide  enough 

information  to  get  better  damping  calculations. 

Finally,  it  should  be  noted  that  the  autoregressive  identification  scheme  did 
not  perform  as  well  as  the  ERA  or  ERA/MME  techniques.  This  is  because  the 
observation  windows  were  too  short.  As  the  windows  get  shorter,  the  autore¬ 
gressive  technique  becomes  an  approximation  to  the  ERA  approach  with  only  a 
modest  amount  of  averaging.  In  the  future,  as  longer  data  sets  are  sampled  the 
autoregressive  approach  should  perform  better. 


window 

length 

sec. 

no.  samples 

Day91008 

Day91010 

24.5 

49.0 

310  j 

616 

Table  5.1:  Length  of  Data  samples  for  Observation  window 
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Figure  5.1:  Measured  Relative  velocity  Day91008 
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Figure  5.2:  Measured  Relative  velocity  Day91010 


S3 


Observation  time  in  seconds 


Model 

Mode 

Hz. 

ERA 

Hz. 

std. 

dev. 

Hz. 

ERA/MME 

Hz. 

std. 

dev. 

Hz. 

AR/model 

reduction 

Hz. 

0.1298 

0.3238 

0.1225 

0.3245 

1  0.5219 

0.0005 

0.0003 

0.0025 

0.1226 

0.3248 

0.5201 

0.0005 

0.0006 

0.0030 

0.1226 

0.3243 

0.5220 

— - 1 

Table  5.2: 

:  Day910l 

IS  Identification  Results:  Moda. 

1  Frequencies 

Model 

Mode 

ERA 

% 

std. 

dev. 

ERA/MME 

% 

std. 

dev. 

AR/model 

reduction 

% 

0.1298 

0.3238 

-0.1244 

1.4477 

4.S917 

0.2489 

0.0410 

0.2984 

-0.3518 

1.1137 

4.5929 

0.3027 

0.0473 

0.2631 

-0.04/0 

1.1583 

4.2525 

Table  5.3:  Day9l008  Identification  Results:  Damping  Factors 


S4 


Model 

Mode 

Hz. 

ERA 

Hz. 

std. 

dev. 

Hz. 

ERA/MME 

Hz. 

std. 

dev. 

Hz. 

AR/model 

reduction 

Hz. 

0.0191 

0.1298 

0.3238 

0.0208* 

0.1244 

0.3312 

0.5115 

0.0033+ 

0.0010 

0.0009 

0.0061 

0.0210  * 
0.1245 
0.3320 
0.5120 

0.0020* 

0.0011 

0.0008 

0.0064 

0.1226 

0.3243 

0.5220 

Table  5.4:  Day91010  Identification  Results:  Modal  Frequencies 


Model 

Mode 

ERA 

% 

std. 

dev. 

ERA/MME 

% 

std. 

dev. 

AR/model 

reduction 

% 

0.0191 

0.129S 

0.3238 

1.8385* 

2.3233 

2.1114 

10.4537 

1.8277* 

0.5182 

0.2312 

1.2157 

1.3937* 

2.1029 

2.0058 

10.7973 

2.2607* 

0.3007 

0.2971 

1.0233 

-0.0470 

1.1583 

4.2525 

Table  5.5:  Day91010  Identification  Results:  Damping  Factors 
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5.2.2  The  Linearity  Assumption 


The  identification  algorithms  are  based  on  a  linear  model  assumption.  While 
it  can  be  assumed  the  minimum  model  estimation  optimal  input  compensates 
for  small  nonlinear  components,  the  extension  to  nonlinear  cases  is  weak.  Thus, 
the  assumption  of  “near”' linear  behavior  is  important.  A  qualitative  study  of 
this  is  done  by  integrating  the  observed  response,  and  then  calculating  the  peak 
displacement  on  the  spacecraft  body. 

The  detected  motion  was  filtered  to  isolate  each  modal  component.  This 
filtered  signal  was  integrated  to  calculate  peak  boom  tip  displacement.  The 
peak  tip  displacement  was  then  divided  by  the  modal  frequency  displacement 
ratios  in  Table  3.2.  These  approximate  peak  system  displacements  are  shown 
in  Tables  5.6  and  5.7.  This  table  shows  that  the  peak  system  displacements 
should  be  small  and  the  linear  assumption  is  valid.  Additionally,  considering 
the  forcing  function  used  (the  deployment  and  retraction  of  the  boom),  the 
vibrations  should  be  linear.  (Note:  No  information  on  system  displacement  at 
.52  Hz  exists). 


mode 

hz 

peak  boom  tip 
Displacement  (mm) 

Displacement 

Ratio 

peak  system 
Displacement  (mm) 

.1298 

0.9 

0.0063 

«  150 

.3238 

1.7 

0.0565 

«  31 

.5200 

0.9 

Table  5.6:  Peak  System  Vibration  Displacements:  Day91008 
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mode 

hz 

peak  boom  tip 
Displacement  (mm) 

Displacement 

.Ratio 

peak  system 
Displacement  (mm) 

.0190 

.9 

0.0242 

«  36  | 

.1298 

1.0 

0.0063 

«  164 

.3238 

1.8 

0.0565 

«  32 

.5200 

.5 

Table  5.7:  Peak  System  Vibration  Displacements:  Day91010 


5.3  Convergence  of  the  Hankel  operator 

As  mentioned  in  the  previous  chapter,  the  identification  algorithm  assume  a 
linear  impulse  response  is  provided.  Given  such  a  response  and  assuming  nucle- 
arity  (see  previous  chapter),  the  Hankel  operator  constructed  from  the  impulse 
response  will  have  a  bounded  nuclear  norm,  Hankel  and  Frobenius  norms.  How¬ 
ever,  the  condition  of  added  noise  may  nullify  these  conditions.  Thus,  a  test  of 
these  norms  on  the  raw  data  for  both  tests  was  performed.  A  data  sample  rate 
of  3.1  Hz  was  used  for  each  filtered  data  set  on  day  91010  and  a  data  sample 
rate  of  6.2  Hz  was  used  for  day  9100S.  Figures  5.6  through  5.17  plot  the 
analysis  results  for  each  data  set.  It  is  clear  that  these  norms  do  not  appear  to 
be  converging.  However,  is  must  be  noted  that  both  Day91008  and  Day91010 
have  limited  sample  periods  and  the  conclusions  drawn  are  based  on  limited 

information. 

1 .  Figures  5.6  and  5.12  show  the  nuclear  norm  for  progressively  larger  Hankel 
matrix  operators. 
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2.  Figures  5.7  and  5.13  show  the  Cauchy  sequence  of  nuclear  norms  for  pro¬ 
gressively  larger  Hankel  matrix  operators.  The  sequence  equals  ||Fn+i  ||jv  — 

||rn|U. 

3.  Figures  5.8  and  5.14  show  the  Frobenius  norm  for  progressively  larger 
Hankel  matrix  operators. 

4.  Figures  5.9  and  5.15  show  an  upper  bound  on  the  Cauchy  sequence  of 
Frobenius  norms.  The  upper  bound  on  the  sequence  is  ||rn+i||f  -  ||rn||jr 
(Triangle  Inequality). 

5.  Figures  5.10  and  5.16  show  the  Hankel  norm  for  progressively  larger 
Hankel  matrix  operators. 

6.  Figures  5.11  and  5.17  show  an  upper  bound  on  the  Cauchy  sequence  of 
Hankel  norms.  The  upper  bound  on  the  sequence  is  ||rn+1||tf  -  ||rn||tf 
(Triangle  Inequality). 

These  plots  show  that  in  spite  of  the  noise  elimination  provided  by  taking  the 
mean  value  over  time  slots  and  filtering,  the  constructed  Hankel  operator  is  not 
converging  to  a  sufficiently  small  bound  within  the  time  period  measured.  This 
could  be  due  to  either  the  continued  presence  of  noise  or  nonlinear  behavior 
of  the  system.  The  latter  of  these  two  sources  is  more  difficult  to  deal  with, 
however  as  discussed  in  the  previous  section  the  linear  assumption  should  be 
reasonable. 
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In  figures  5.18  through  5.29  the  same  analysis  was  done  on  the  minimum 
model  error  data  (see  previous  chapter  for  MME  description).  As  discussed  in 
the  previous  chapter,  this  data  is  constructed  optimally  from  a  linear  model 
such  that  the  exogenous  input  is  minimized.  In  this  case  the  exogenous  input 
can  be  considered  to  be  nonlinear  components  entering  into  the  system.  These 
plots  show  that  the  nuclear  and  Frobenius  norms  are  beginning  to  converge.  By 
allowing  the  target  covariance  used  in  the  optimization  to  increase,  the  norms 
will  converge  even  faster.  The  target  covariance  used  was  2.4  mm/sec  as  in  the 
previous  section. 
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Figure  5.7:  Cauchy  Sequence  of  Nuclear  Norms  for  Truncated  Hankel  Operator 
Day91008 
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Hankel  Operator  Day9100S 
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Hankel  Operator  Day91008 
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Hankel  Operator  Day91010 
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Figure  5.23:  Upper  Bound  of  Cauchy  Sequence  for  Hankel  Norm  of  Truncated 
MME  Hankel  Operator  Day91008 
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Figure  5.26:  Frobenius  Norm  Truncated  MME  Hankel  Operator  Day91010 


Figure  5.27:  Upper  Bound  of  Cauchy  Sequence  for  Frobenius  Norm  of  Truncated 
MME  Hankel  Operator  Day91010 
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5.4  Conclusions 


Considering  the  brevity  of  the  observation  windows,  the  experiment  was  ex¬ 
tremely  successful.  The  experiment  accomplished  the  following  goals: 

1.  Structural  modes  were  verified  independently  using  the  actual  spacecraft 
in  its  true  environment.  Three  of  the  first  four  modes  were  detected. 
These  results  are  impressive  compared  to  previous  experiments,  specifi¬ 
cally,  the  19S5  Solar  Array  Flight  Experiment  (SAFE).  The  SAFE  exper¬ 
iment,  which  was  a  two  day  experiment  done  on  board  the  space  shuttle, 
tested  the  modal  properties  of  a  deployed  truss  with  an  attached  solar 
panel.  The  SAFE  experiment  only  detected  three  modes.  Comparing  the 
overall  experimental  cost,  this  experiment  done  on  the  LACE  spacecraft 
was  an  exceptional  success. 

2.  Algorithmic  approaches  were  tested  with  an  actual  experiment.  Algo¬ 
rithms  provided  robust  identification.  The  identified  modes  were  tightly 
grouped  and  agreed  with  predicted  modes.  A  variety  of  identification  tech¬ 
niques  were  compared  on  simulated  data.  The  Hankel  operator  technique 
performed  better  than  any  other  approach,  given  the  type  of  signal  and 
the  short  length  of  observation. 

3.  A  new  laser  application  was  verified.  This  experiment  showed  that  laser 
tracking  could  be  accomplished  and  provide  data  on  which  analysis  could 
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be  performed.  This  approach  offers  alternative  experimental  approaches 
to  structural  identification  for  low  orbit  satellites. 

4.  A  system  engineering  approach  was  applied  to  a  complex  problem,  whose 
technological  solution  spanned  many  disciplines.  A  requirement  analysis 
review  showed  the  need  for  large  space  systems.  Given  this  requirement, 
the  structural  identification  is  a  necessary  component  of  on-station  health 
monitoring  and  real  time  control  algorithms.  •  . 

This  last  point  is  very  important.  The  problem  of  spacecraft  dynamics  is 
generated  by  increased  demands  from  power  considerations  and  the  increased 
size  of  space  structures.  This  experiment  was  designed  to  show  that  low  cost  low- 
complexity  on-board  hardware  combined  with  algorithm  analysis  can  satisfy  the 
identification  process  and  aid  in  the  progression  toward  real-time  identification 
and  control. 
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